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CHAPTER 0 : INTRODUCTION 


Segmentation is a major technique used in the quantitation in magnetic 
resonance imaging (MRI). This chapter recalls the basics of MRI, and contains a brief 
survey of the literature on the segmentation of images. This thesis is a study of 
segmentation of MRI images using different filters, such as Gaussian and Butterworth 
filters and discrete convolutions with Korovkin, Jackson, generalized Jackson, and De 
La Vallee-Poussin kernels, followed by maximum likelihood estimation (MLE) based 
segmentation algorithms developed both for single contrast, as well as multisequence 
images. 

Magnetic resonance imaging is a noninvasive diagnostic technique that 
produces images of internal body tissues, which is based on nuclear magnetic 
resonance of atoms within the body induced by the application of radio waves. MRI 
dominates over other medical imaging modalities due to its excellent soft-tissue 
contrast, high spatial resolutions and non-invasiveness. MR images are degraded due 
to, among others, the inhomogeneities in the radio-frequency field (known as bias 
field intensity) and the presence of multiple tissues within a pixel that leads to a 
partial volume effect. Presence of many neurological diseases alters the shape and the 
volume of brain parenchyma and cerebrospinal fluid (CSF) regions. So, their 
segmentation, and the resulting classification of tissues and pathologies are important 
in many medical imaging applications. Segmentation helps in the classification and 
quantitation of normal tissues, as well as different pathologies. Segmentation also 
helps in obtaining computerized anatomical structures, and in guiding surgery. 

An image segmentation is defined as an exhaustive partitioning of an input 
image into regions, each of which is considered to be homogeneous with respect to 
some image property of interest (e.g., intensity, color or texture). Segmentation is one 
of the major steps in the analysis of processed image data. It has a different 
significance in different areas. Like its main aim in aerial photography is to divide a 
photographic or geographical or satellite image into parts that have a strong 



correlation with objects or areas of the real world contained in the image. In medical 
images, like ultrasound, computerized tomography, magnetic resonance images, the 
main aim of segmentation is to separate the region of interest from rest of the image. 
Segmentation of magnetic resonance images is complex due to various kinds of 
degradations, and as such no single method seems to be applicable to all kind of 
images. 

In MRI, Segmentation helps in the classification of normal tissues along with 
the pathologies present within them. It also helps in the quantitative analysis that 
provides a better knowledge of the growth and diagnosis of pathology during an 
interval of time and this information could be used in further treatment. It is useful in 
the understanding of physiology of brain and other parts of the body by enabling a 
study of the various chemical and biological environments of the tissues. 
Misclassification of pixels into different tissue classes due to the presence of noise 
could give improper segmentation which consequently would lead to poor tissue 
quantitation. The proper choice of filter to reduce noise would help in minimizing 
these misclassifications and the corresponding errors. 


0.1. BASICS OF MAGNETIC RESONANCE IMAGING 

Precessional frequency: Magnetic nuclei, when placed in an external 
magnetic field, Bo, possess a unique precessional motion with Larmor frequency, ft, 
which is directly proportional to both the strength of the external magnetic field and 
gyromagnetic ratio, y , of the magnetic nuclei and is given by 

The gyromagnetic ratio is defined as the ratio of the magnitude of the magnetic 
moment of the nucleus and the magnitude of its spin angular momentum. 

Resultant Magnetization Vector; Resultant magnetization vector M of 
nuclei is the net sum of all the magnetic moment vectors which align in the direction 
of external magnetic field. Bo, on the application of Bo. 



Radio-frequency Magnetic Field: The perpendicular radio-frequency field is 
produced using a brief application of an alternating electric current through a coil. 

Whenever radio-frequency magnetic field Bi, perpendicular to M, is imposed, 
it results in a precession of M around Bi tipping away its orientation from its original 
position, due to which M acquires a component perpendicular to that of Bo causing it 
to process around Bo also. M just wobbles about its original direction as long as Bi 
continues to change at an arbitrary frequency. The angle of tip continues to become 
large in resonance, which is caused when the frequency of Bi equals the Larmor 
frequency due to the field Bo- Besides its longitudinal component Ml, the transverse 
component Mj of M is produced as M is tipped from its original orientation (Figure 
0.1) following radio-frequency pulse (C.L. Partain, R.R. Price, J.A. Patton, M.V. 
Kulkarni and A.E. James Jr. [1988]). . 


z 

Bo i 



Figure 0.1: The Larmor precession geometry M traces out a circle parallel to xy plane and its 
component Mt traces out a similar circle in the xy plane, while other component Ml remains fixed. 


A sustained pulling of M away from zero can occur only if Bi is always 
perpendicular to M which makes Bi to rotate in the same direction as M at the rate fi 
(Figure 0.2). With Bi rotating at this rate, M remains perpendicular to it while Larmor 
processing around it at the rate; 



/./ = 


2n 


B 


1 ’ 


where /sf is the spin flip precession rate. 



Figure 0.2: Left: Bi is stationary along the x-axis. It pulls M clockwise away from z-axis, but Bo also 
pushes M toward the xz plane. Mt has rotated through the angle changing the angle between Bi 
and M. Right. Bi has been rotated as the same rate as Mt: this keeps the angle between M and Bi 
fixed, and hence the force on M is the same, pulling M away from Bo even further. 


During relaxation, Bi is absent, Ml grows from its spin-flipped value toward 
an equilibrium value of M in a fully relaxed sample. The rate of relaxation is 
proportional to the difference between Ml's current value and its equilibrium value at 
full relaxation. The time constant of this relaxation is called Ti, and the growth rate of 
Ml is referred to as Ti relaxation, also termed as spin-lattice relaxation. At the same 
time, Mt decreases to zero as it rotates at the rate fi around xy plane. The rate of 
decrease of Mj to its equilibrium value of zero is proportional to the current 
magnitude of Mt. The time constant of this relaxation is called T 2 , and the decay rate 
of Mt is referred to as T 2 relaxation (also termed as spin-spin relaxation) (C.L. 
Partain, R.R. Price, J.A. Patton, M.V. Kulkarni and A.E. James Jr. [1988]). 



The maximum signal in the receiver coil, following 90° pulse starts decreasing 
owing to relaxation. This property of signal is known as the free induction decay 
(FID). Pulse cycle is' a repeating unit, which is composed of a series of one or more 
radio-frequency pulses with a measurement of one or more MR signals. Pulse 
sequence is a series of pulse cycles. TR (Time to Repeat) is the time interval between 
two successive pulse cycles, measured in milliseconds. TE (Time to Echo) is the time 
from one pulse to the measurement of MR signal, which is also measured in 
milliseconds. Ti weighted images owe their contrast mainly to Ti relaxation 
properties of the tissues. They use a short TE, and a short TR. T 2 weighted images, 
with contrast predominantly due to the Ta relaxation properties of the tissues, are 
produced by using long TE and a long TR. Proton Density weighted image is 
produced using long TR, with short TE (A.L. Horowitz [1989]). Their contrast 
depends on the concentration of water molecules in the tissues. 

Magnetization transfer (MT) contrast imaging is used to suppress background 
tissues to enhance vessels and certain disease processes. In this process, prior to the 
conventional pulse sequence procedure, a pre-saturation pulse is used to saturate the 
bound protons and promote an exchange of some of this saturated magnetization onto 
the free protons. This results in reduced signal from the free protons, e.g., gray and 
white matter lose 30-40% of their signal when an MT pulse sequence is utilized. 

The spin echo (SE) sequence can be written as (90° - TE/2 - 180° - TE/2 - 
measure - TD)n (TD is a delay time). Since initially M is aligned with Bo, Ml = M 
and Mt = 0. Following the 90° pulse. Ml = 0 and Mt = M, with both Ti and T 2 
relaxation beginning and allowed to continue (as FID) for a time TE/2. The relaxation 
equations are given as: 

Mi (0 = M (1 - (0)e""^' 

Mj (t) = Mj (True Tj decay) 

MT(t)=MTi0)e~‘'^^' (FID decay). 

Here Ml( 0) is the value of Ml following a radio-frequency pulse. Ml is completely 
relaxed at the start of the first sequence, but only partially relaxed at the end of the 
sequence, which is a recovery from the 90° and 180° pulses. Ml ends with same 



partially relaxed value on second and subsequent sequences with which it started after 
the first sequence and the signal intensity can be written as: 

0.2. SEGMENTATION 

In magnetic resonance imaging, the segmentation methods are broadly 
classified into two categories based on a single image contrast and multi-sequence 
images (L.P. Clarke, R.P. Velthuizen, S. Phuphanich, J.D. Schellenberg, J.A. 
Arrington and M. Silbiger [1993], L.P. Clarke, R.P. Velthuizen, M.A. Camacho, J.J. 
Heine, M. Vaidyanathan, L.O. Hall, R.W. Thatcher and M.L. Silbiger [1995]). 
Methods corresponding to a single image are further classified into two categories - 
edge based, and region based methods. Edge based methods consist of edge detection, 
boundary tracing and thresholding; whereas region based methods are seed growing, 
merge and split algorithms and template matching. For multi-sequence images, 
supervised and unsupervised algorithms are mentioned in the literature. 

Edges in camera photographs are due to a change in some physical properties 
like surface reflectance, the geometry and the intensity of an object (A.K. Jain 
[1995]). The well- known edge detector, Sobel operator, is defined as the square root 
of the sum of squares of gradients in x- and y-directions. Edges are also detected by 
finding zero-crossings in the convolutions of the image with Laplacian-of-Gaussian 
masks (A. Huertas and G. Medioni [1986]), or by Marr-Hildreth operator (R.M. 
Haralick [1984], E.C. Hildreth [1983]) which can be approximated by Difference-of- 
Gaussians that seems to be good to obtain the edges for large regions. Marr-Hildreth 
operator is used along with morphological filtering (M. Bomans, K-H. Hohne, U. 
Tiede and M. Riemer [1990]), which requires manual labeling and editing of the 
regions to generate satisfactory 3-D displays. The edge detection algorithms generally 
suffer from incorrect detection of edges due to noise, over segmentation due to the 
presence of inhomogeneity within the region and under segmentation due to 
inaccurate boundaries between the adjacent regions. 



Boundaries are obtained by joining the edge element at a certain pixel to that 
present in its first order neighborhood or in second order neighborhood. First order 
neighborhood consists of four horizontal and vertical pixels, while second order 
neighborhood has first order neighborhood along with four diagonal pixels. In contour 
following algorithms (G.P. Ashkar and J.W. Modestino [1978]), the pixel is viewed as 
having a square-shaped boundary, and the boundary is traced by following the edge 
pixels using graph searching methods, where a pixel acts as a node in a graph (U. 
Montanari [1971]). In some cases, boundaries are obtained using Fourier descriptors 
(A. Chakraborty, L.H. Staib and J.S. Duncan [1996], L.H. Staib and J.S. Duncan 
[1992]) and B-splines (A.K. Jain [1995]). 

Thresholding produces a segmentation, which gives all the pixels, either 
belonging to objects of an interest or background of an image. This is mainly based on 
the shape of the histogram of image data, such as number of non-overlapping peaks 
present in histogram plot according to which single or multiple thresholds can be 
predefined and the thresholding is known as a simple thresholding or multiple 
thresholding (J.C. Russ [1992]). Optimal thresholding is obtained by the 
approximation of the histogram of an image using a weighted-sum of two or more 
probability densities with normal distributions. These algorithms in magnetic 
resonance imaging are also affected by inhomogeneity present within region of 
interest. 

Seed growing is a region-based method, which selects the seed and predefined 
threshold, depending on the region of interest. It then examines the pixels surrounding 
the seed on the basis of intensity or some textural properties of seed and includes the 
neighboring pixel into the region if it comes within a pre-decided threshold (R.C. 
Gonzalez and R.E. Woods [1999]). The inclusion of pixels is mainly according to first 
and second order neighborhood. Each new added pixel then becomes a seed and the 
process continues till there is no inclusion of new pixel into the region. This method is 
helpful in segmenting binary images and it is mainly used as a post processing 
technique to operate on segmented image obtained from some other methods. 



Region splitting and merging (R.C. Gonzalez and R.E. Woods [1999]) is an 
iterative procedure which starts with the splitting of an image into disjoint sets and 
then the merging takes place and stops when no further splitting and merging is 
possible. Splitting of a region R on an image into its n subregions Ri, Rz, ... , i?n, and 
merging between any two regions are done based on the following conditions: 

(1) U-=iR/=/?, 

(2) Ri is a connected region, i = 1, 2, ..., n, 

(3) Ri nRy= <j> for all i and j, i ^ j, 

(4) P(Ri) = TRUE for j = 1, 2, ..., n, 

(5) P{Ri u Rj) = FALSE for i ^ j, 

Here P(Ri) is a logical predicate over the points in set R,- and 0 is the null set. 

Template matching is based on the matching of a small section, either of 
certain tissue or abnormality, with the given image (A.K. Jain [1995]). This generally 
gives location of abnormality and fails to provide any kind of distinction either in 
terms of boundary or in terms of regions between pathology and normal tissues. 

Supervised methods require initial training data to begin with, which is mainly 
based on rough segmentation or prior knowledge of some features of an image. One 
of the famous non-parametric supervised classifiers, k-Nearest neighbor (M. 
Vaidyanathan, L.P. Clarke, C. Heidtman, R.P. Velthuizen and L.O. Hall [1997], L.P 
Clarke, R.P. Velthuizen, S. Phuphanich, J.D. Schellenberg, J.A. Arrington and M. 
Silbiger [1993]), uses a sample set of pixel vectors Z and assigns them to different 
tissue classes. The unlabeled training pixel vectors X are assigned according to the 
Euclidean distance between X and the labeled pixel vectors Z. The posteriori 
probability for X is then estimated from the frequency of the labels k-NNs and the 
tissue class label is assigned to X based on the maximum posteriori probability. Fuzzy 
c-means (D.L. Pham and J.L. Prince [1999]) is a well-known supervised algorithm for 
segmentation. 

Unsupervised methods automatically find the structure or pattern present 
within the image and one such well known method is clustering. It is a region based 



method which uses the spatial co-ordinates of the centroids of features which are 
usually gray level intensity, texture, etc., and sorts through them to locate the nearest 
neighbor features for each feature present. It is based on center-to-center distances 
and then constructs a distribution image of the intensity of nearest neighbor distances. 
Some unsupervised techniques could be further processed via multi-spectral analysis 
along with region analysis (M.C. Clark, L.O. Hall, D.B. Goldgof, R. Velthuizen, F.R. 
Murtagh and M.S. Silbiger [1998]) to get the required segmentation. 

Some algorithms (M.C. Clark, L.O. Hall, D.B. Goldgof, R. Velthuizen, F.R. 
Murtagh and M.S. Sibiger [1998]) integrate knowledge-based techniques with multi- 
spectral analysis. Initial segmentation is obtained by an unsupervised clustering 
algorithm following the multi-spectral histogram analysis and classification of 
segments based on region analysis. Knowledge of tissue intensity properties and 
intensity inhomogeneities are used to correct and segment MR images via 
expectation-maximization algorithm (W.M. Wells III, W.E.L. Grimson, R. Kikinis 
and F.A Jolesz [1996]). 


Segmentations based on normal mixture models of different kinds are also 
mentioned in the literature. The pixel intensity either follows univariate normal 
distribution in case of a single image or multivariate normal distributions in case of 
multiple images. 


In Y. Zhang, M. Brady and S. Smith [2001], it is assumed that the pixel 
intensity follows a Gaussian distribution with parameters 0^ given the 

class label x,- =1, 


P(yi\Xi)= ^exp.- _ 2 

(J I y 271 2(7 ^ 

and the joint likelihood probability, based on the conditional independence 
assumption of y is given by 


n. , s — . I N 1 (yi-Pxf , 

P(y|x) = _^p(y,lx,)= _^/^expX - - a -logcr,^ 

(27r) /e5^ 2cr^. ^ 




ieS 



where 5 is a two dimensional rectangular lattice and each pixel is characterized by an 
intensity value y- and the labeling of S is denoted by x, where x,-, {ie S) is the 

corresponding class label of pixel i. The random field X = {Z,-, is 5}is an 
underlying Markov random filed (MRF) assuming values in a finite state space L with 
probability distribution P(x) = exp(-f/(x)) (J. Besag [1974]), where 

Z = ^ normalizing constant called the partition, and U(x) is an 

energy function of the form t/(x) = which is a sum of clique potentials 

(x) over all possible cliques C. A clique C is defined as a subset of sites in S in 

which every pair of distinct sites are neighbors, except for single-site cliques (S. 
Geman and D. Geman [1984]). Labeling of class in an image is done according to the 
maximum a posteriori (MAP) criterion, 

X = arg max(P(y/x)P(x) } . 

xeX 


A. Lundervold and G. Storvik [1995] propose a method for segmentation of 
brain parenchyma and cerebrospinal fluid, where head and brain are divided into four 
major regions and seven different tissue types. Each region is associated with a finite 
mixture of its constituent tissue types and each tissue is modeled by a multivariate 
Gaussian distribution. In this paper, initial estimates of tissue parameters, means and 
covariances of mixture are obtained by using k-means clustering algorithm. The 
boundary of each region is extracted using contour detection algorithm followed by 
the classification and distribution of tissues in each of the regions. This segmentation 
method is restricted to images where the brain parenchyma and CSF spaces form 
connected regions. 

The process of segmentation is to find x which represents the correct tissue 
class at each voxel site given by image y (J.C. Rajapakse, J.N. Giedd and J.L. 
Rapoport [1997]). To find the maximum a posteriori (MAP) estimation from the 
image data, that is, ifx = x* represents the optimal segmentation 

X* = arg max p{x \ y) , 

xgQ 

where Q is the set of all possible segmentations and p(x\y) is the posterior density of 
the segmentation x given the image y. From Bayes' theorem, p(x\y) p{x,y) = p(y[x). 



p{x) due to the independence of prior probability of image p{y), where p{y\>c) is the 
conditional density of a pixel intensity y given segmentation x as: 

'' 2 '' 

, I s 1 

p(y\x) = I — exp^-- 

where Rk denotes the region or the set of all voxel sites belonging to tissue class k, p.k,i 
and Ok^i represent the mean intensity and standard intensity of the class k at pixel site i. 

Comparisons of different segmentation algorithms according to the 
performance on intensity inhomogeneities within each tissue region, along with their 
stability and applications (M. Vaidyanathan, L.P. Clarke, C. Heidtman, R.P. 
Velthuizen, and L.O. Hall [1997], L.P. Clarke, R.P. Velthuizen, S. Phuphanich, J.D. 
Schellenberg, J.A. Arrington, and M. Silbiger [1993], D.L. Pham, and J.L. Prince 
[1999]), like maximum likelihood method, clustering, artificial neural networks (S.C. 
Amartur, D. Piraino and Y. Takefuji [1992], W.E. Reddick, J.O. Glass, E.N. Cook, 
T.D. Elkin and R.J. Deaton [1997]) are also reported in the literature. 



0.3. TEST IMAGE SETS 


All numerical simulations in this thesis are done on the data obtained from 
Department of Radio-diagnosis, Sanjay Gandhi Postgraduate Institute of Medical 
Sciences, Lucknow. Magnetic Resonance Imaging of patients are performed on a 1.5 
T super-conducting system using circularly polarized head coil. Conventional Proton 
density (PD), T 2 (TR/rE/l,2/n = 2200/12, 80/1) and Ti (1012/14/2) weighted SE 
images are acquired in axial plane using 256x256 matrix size, 0.1mm inter slice gap 
and 5mm slice thickness. We have considered three sets of images for our numerical 
simulations in subsequent chapters. Figure 0.3 shows set I which consists of T 2 
weighted image and Figure 0.4 shows set II which consists of Proton density, T 2 and 
Ti weighted images. Set III consists of Proton density, T 2 , Ti and MT SE Ti weighted 
images shown in Figure 0.5. Magnetization transfer Ti weighted imaging is 
performed with identical parameters as for Ti weighted image except for an off- 
resonance pulse. Sets I and II are the images of different cross-sections of a patient 
with edema located in the region of brain parenchyma. Set III is a case of post- 
traumatic epilepsy with perilesional gliosis. 



T 2 weighted Image 


Figure 0.3: Test Image Set I 




T 2 weighted Image 


Ti weighted Image 


Figure 0.4; Test Image Set II 







Proton Density weighted Image T 2 weighted Image 



T 1 weighted Image MT SE Ti weighted Image 


Figure 0.5: Test Image Set III 







0.4. CONTENTS OF THE THESIS 

A detailed study of the degree of approximation properties of the Marr- 
Hildreth operator and other related convolutions with the Gaussian function and the 
resulting segmentations has been done in Chapter I. 

Asymptotic estimates for Marr-Hildreth operator are deduced in section 1.3 
which are as follows; 

Theorem 1.3.1. 

If I(X) is 4-differentiable at a point X, 

lim^^o(T-^{l(X)0V^G(X,<y)-V^I(X)}=V^I(X)/2, 

and this relation holds uniformly in X 6 5 if I(X) is uniformly 4-differentiable in S. 
If 1(X) is 2p-differentiable at a point X, 

X 0 /(Z ) / dx^‘ dy ) + 0(CT , 

where the sums run over 2 < k < p, I, m, n > 0, with l+m+n = k. Moreover, the 
relation holds uniformly in X e S if /(X) is uniformly 2p-differentiable in S. 

Theorem 1.3.2. 

If /(X) is 2p-differentiable at a point X, 

2cT-"[/(X)®G(X,cr)-/(X)]-V2/(X) = 25 ^^a'^'^-'^ 2 ;,,,„(( 2 ^)-( 2 m)!!( 2 n)!!)-‘ 

X O^^/fX) / dx^‘dy^'”dz ^'' ) + ) , 

where the sums run over 2 < k < p, I, m, n > 0, with l+m+n = k. This relation holds 
uniformly in X £ if /(X) is uniformly 2p-differentiable in S. In particular, if /(X) 
is 4-differentiable at a point X, 

lim^_o(T-2[2cr-2[/(X)®G(X,cT)-/(X)]-V2/(X)] = V"/(X)/4, 
and the relation holds uniformly in X e S if /(X) is uniformly 4-differentiable in S. 


Section 1.4 deals with the computation of gradients of arbitrary order. 



Theorem 1.4.2. 

Let /(Z)be a bounded continuous function, (^-differentiable at a point X. If Q(u,v,w) 
is any polynomial of total degree q in u, v, w, then with QiD) denoting the 
differential operator Q(d /dx,d/ dy, d / dz), 

lim^^o /(X)0[(2(£>)G(X,(7)] = Q(D)I(X), 
the relation holding uniformly in X 6 S' if /(X) is uniformly ^-differentiable in S. 
Moreover, if /(X) is 2p+^-differentiable atX, then 

Q(D)[I(X)® G(X,ct) - /(X)] = 'Zx,My ((2A)!!(2At)!!(2v)!!)-' 

x{d^'‘[Q(D)l(X)]/dx^^dy^^dz^'')+o(a^P), 
where the sums run over 2 < k<p, X,/u,v>0, with A + ^ + v = ^ , and the relation 
holds uniformly in X 6 S if /(X) is uniformly 2p-i-^-differentiable in S. 

In Section 1.5, the degree of approximation of C^-functions by a Gaussian 
convolution is studied. The following are the main results: 

Theorem 1.5.4. 

Let /(X)€ C, and y/e S. Then: 

(i) /(X)0 G(X,<t)-/(X) =O(i/r(o-^)),cr^0,iff X(/;r) = 0(vr(r^)), r-»0; 

(ii) /(X)0G(X,cr)-/(X) =o(v(c7^)),cr-^0,iff X(/;r) = o(i/r(r^)),r->0; and, 

(iii) /(X)0G(X,cr)-/(X) ~ X(/;(t), if X(/;r'^^)6 S . 

Here, and in the sequel, S denotes the class of continuous non-decreasing functions 
y/(t) satisfying 0 = y/iO) < y/(t),t > 0 such that limj,^o sup o<, < 1/2 = 0 • 

Theorem 1.5.6. 

Let Q(D)1(X)eC, and y/eS. Then: 

(i) I(X)®[Q(D)G(X,(y)]-Q(D)I{X) = Oiy/ a -^0, 
iff K(QiD)I(X);t) = O(y/it^)),t^0-, 

(ii) /(X)0[(2(r>)G(X,CT)]-!2(Z))/(X) =o(y/{(j^)),(j-^0. 



iff K{Q{D)l{X)-,t) = o{y/it^)),t-^ 0 ; and. 

(iii) I(X)®[Q{D)G(X,<j)]-Q(D)IiX) ~K{r,a), if K{Q{D)l{Xy,t^'^)e S. 
Here Q{u,v,w) is a polynomial in u, v, w, as before, and Q(D)=Q(d/dx, d/dy, d/dz)- 

Saturation of convolutions with Gaussian filters has been discussed in section 
1.6. Let the radial average Ip (X) of a function l(X)e C be defined as; 

lp{X) = (47rp2)-‘ “ l{T)dA{T), p > 0, 

^ ••ACa.p) 

where A{X,p) denotes the surface of the sphere about the point X of radius p, T is a 
general point and dA{T) stands for the area element on the surface about the point T. 
We call 1{X) harmonic if lp(X) = 1(X), Xe R^, p > 0. In particular, the following 
results are obtained: 

Theorem 1.6.1. 

Let I{X)e C. Then: 

(i) /(X)0G(X,cr)-/(X)| = O(a^),(7-^O, iff /,(X)-/(X)| =O(r^),r^0; 
and, 

(ii) \liX)®G{X,(j)-liX) =o((T^),cr -^0, iff /(X) is harmonic. 

Corollary 1.6.2. 

Let (2(Z))/(X)e C. Then: 

(i) /(X) 0 [Q(D)GiX,CT)] - QiD)IiX)\\ = 0{a\ -> 0, 
iff Q(D)I,(X)-Q(D)1(X) =O(r"),r^0;and, 

(ii) IiX)<S>[Q(D)G(X,(7)]-QiD)IiX) = a ^0, 

iff Q(D)I(X) is harmonic. 

In section 1.7, the following aysmptotic error estimates are obtained for 2-D 


Gaussian filter: 



Theorem 1.7.1. 

If I(X) is 4-differentiabIe at a point X, 

lim^^oCr-^{l(X)®V^G(X,a)-V-l(X)}=V^I(X)/2, 

and this relation holds uniformly in X e 5 c R^, if 7(Z) is uniformly 4-differentiable 
in 5. If I(X) is 2p-differentiable at a point X, 

l(X)<S>V^G(X,a)-V^IiX) = J^^(2k)a^^^-^^Y^^J(2l)l\(2m)\\)-^ 

x(d^‘‘l(X)/dx^‘dy ) + o(cr , 

where the sums run over 2 < k < p, I, m> 0, with l+m = k. Moreover, the relation 
holds uniformly in X 6 5 if I(X) is uniformly 2p-differentiable in S. 

Theorem 1.7-2. 

If 1(X) is 2p-differentiable at a point X, 

[I(X)®G(X,ct)-1(X)] = (<T^ /2)V^1(X) 

the sums being over 2 < k < p, with I, m > 0, and l+m = k. This relation holds 
uniformly in X € 5 c if/(X) is uniformly 2p-differentiable in S. In particular, if 
I{X) is 4-differentiable at a point X, 

\im^^o^-^^a-^U(X)(8)G(X,cj)-liX)]-V^IiX)}=V^l(X)/4, 
and it holds uniformly in X € 5, if /(X) is uniformly 4-differentiable in S. 

Corollary 1.7.3. 

For the difference of Gaussians (DOG) in two variables: 

/(X)®[G(X,a,)-G(X,CT,.)]-ct/(l-r")V2/(X)/2 

= ct/(1 - /)V''/(X) /8 + o(<t/), cr, ^ 0. 

The two-dimensional gradients of arbitrary order have been considered in 
section 1.8. The following is a pointwise simultaneous approximation property of an 
arbitrary order for the 2-D convolution. 



Theorem 1.8.1. 

Let /(X) is a bounded continuous function, (j-differentiable at a point X, 
lim^^o nX)®[Q(D)G(X,a)] = Q(D)1(X), 
the relation holding uniformly in X e 5 c if /(X) is uniformly ^-differentiable in 
S. Moreover, if /(X) is 2p+^-differentiable atX, then 

2(£>)[/(X) 0 G(X, (T) - /(X)] = X, EA./(2A)!!(2At)!!)-^ 

x(d^^[QiD)l(X)ydx^%^^)+oi(j^n, 

where the sums run over 2< k<p, A,/i>0, with X + p = k , and the relation holds 
uniformly in X g c R^, if /(X) is uniformly 2p+^-differentiable in S. 

Error for smooth gradient functions in plane has been obtained in section 1.9 
and the main results are as follows; 

Theorem 1.9.4. 

Let /(X)g C, and y/e S. Then: 

(i) l(X)®G(X,cr)-I(X) =O(\i/(a^)),<j-^0, iff X(/;r) = O(i/^(r^)),r-^0; 

(ii) 7(X)0G(X,cr)-/(X) = o(i/r(cr2)),a 0, iff X(/;r) = o(i/^(r^)),r-^0; and, 

(iii) /(X)0 G(X,c7)-/(X) ~X(/;(t), if X(/;r'^^)G5. 

Theorem 1.9.6. 

Let Q(D)I{X)e C, and y/eS. Then; 

(i) I(X)®[QiD)G(X,(7)]-Q(D)nX) =O(wi(7^)),(y-^0, 

iff K(Q(D)I(Xy,t) = O(y/(t^)),t^0-, 

(ii) l(X)®[QiD)GiX,(j)]-Q(D)I(X) =o{y/icj^)),(j^0, 

iff KiQ{D)I{X)\t) = o{yf{t‘^)\t 0 ; and, 

(iii) /(X)0[G(D)G(X,a)]-(2(T>)/(X) ~ X(/;c7), if X((2(i9)/(X);r'^")G 5, 

where Q{u, v) , as before, a polynomial in u, v, and Q{D) = Q{d /dx,dl dy) . 



Numerical simulations and discussions for the Gaussian filter have been done 
in section 1.10. Segmentation results with Gaussian filter corresponding to cases a'* = 
64 and 32 are practically useful as the region of pathology is well segmented in these 
cases. Sobel maps for Gaussian filtered images give better description of the edges as 
compared to that of difference magnitudes of Gaussian filtered images. The contours 
obtained with Marr-Hildreth operator for a = 2.0, 2.576 and 3.0 are also useful. 


Chapter, II is a study of Butterworth 1-D filters and its generalizations. These 
filters are characterized as solutions of certain minimization problems. The main 
results are: 

Theorem 2.2.1. 

For feL 2 „^,m = l,2,3,..., and 6 > 0, the trigonometric polynomial T„{x) e Tn 
minimizing 

Q„=i2n)-^ ' j fix)-T„(xf +en-’^Tj"'\x)^]dx, 


is unique and is given by 

(x) = T„(f-,x) = aQ(f)l2 + i^k (/) cos kx + b^ (/) sin kx)/(^ + d(k/n )^'" ), 

where a^. (/), bj^ (/) are the Fourier coefficients of the function f(x) . The map/-> 
Tn is linear, and, T^- f -^0, 


Here the norm / is given by 

1 -^ + 


^k + 


hr! 


1/2 


Theorem 2.2.3. 

Let /e and yre S. Then: 

(i) There exists a constant M < max(2,0), such that 


/-r„ 

(ii) / - r„ = ), n iff / 6 ; 



(iii) f ~Tn - iff /is a constant; 

(iv) / -r„ = 0(v^(n-2'")),n ^ oc, iff ii:2^(/;r) = 0(1/^/^'")), t 0 ; 

(v) I/-7; I = o(J//(n"^"' )),«-»=«, iff A:2;„(/;0 = o(i//^(t^"')),t->0; and, 

(Vi) / -T„ ~ ii:2„,(/;l/n), if 5. 

HereA-,„(/;() = inf^^^^^_,{l/-s|+,^" «L} 

Section 2.3 studies the two dimensional Butterworth filter. 

Theorem 2.3.1. 

For / £ L2 ;j _2 > = 1> 2, 3, ..., and 0> 0, the trigonometric polynomial r„(x,y) e Tn,2 

minimizing 

is unique and is given by: 

T„ (x, y) = T„ {f\x, y) = (/) 1 + 0(| * i / ^ exp + ky y ) }, 

where (/) are the complex Fourier coefficients of the function f{x,y). The map: / 
— > Tn is linear, and, /, - / ^0, n — > 00. 

The linear space invariant map r„ : f -^T„{f',x, y) has r„ = 1 , and it can be 

evaluated by FFT using multiplicative replacements: ^ /[1 + 0(| it | /n)^'"] in 

the frequency domain. The effect of this filtering for different 6 and m and also the 
degree of approximation of/ (x,y) by the polynomials T^(f;x,y) is considered. 

The main direct and inverse result on the degree of approximation 
of /( Jc, y) by r„(x,y) is: 

Theorem 2.3.3. 

Let /£ L 2 „^ 2 ^, and y/G S. Then: 

(i) There exists a constant M < max(2,0) , such that 



f '^n — 2m, 2^/ jl / /z) — > 0, « — 4 oo ; 

(ii) / ~^n| = 0(n n — > «=, iff / e L2^^22m' 

(iii) |/ “T),! = o(/?~^'"), n -4 oo, iffyis a constant; 

(iv) /-T; =0(i/^(n-2'")),n^oo, iff A: 2 m, 2 (/;^) = O(i/^(r2'")),t^0; 

(v) =o(i//'(n“^'”)),n-4oo, iff ij: 2 ^_ 2(/;0 = o(v/'(r^'”)),t->0;and, 

(Vi) / -T„ ~ ii:2m,2(/;l/«)> if 5. 

Section 2.4 studies the three dimensional Butterworth filter, beginning with 
the following characterization: 

Theorem^2.4.1. 

For / e L 2 ;j _3 , m = 1, 2, 3, ..., and 0> 0, the trigonometric polynomial T„{x,y,z) 6 
Tn ,3 minimizing 

e„ = (2®) { / (^. y. z) - t; u y, z)p 

^ n~”'d'% {x, y, z) / dx‘dy^dz‘ ^ ^dxdydz, 
is unique and is given by: 

(x, y, z) = T„ (f-,x, y, z) = 2i<|q<„ (/) 1 + ^(l ^ I / ' exp + kyy + k^z)}, 
where c^(/) are the complex Fourier coefficients of the function f(x,y,z). The 
map:/— > Tn is linear, and, T^- f ^0, 

As in the 2D case, the linear space invariant map T„: f T„(f',x,y,z) has 
r„ =1, and can be evaluated by FFT using multiplicative replacements: 

l\l + 6(\k\l in the frequency domain. The main direct and inverse 

result on the degree of approximation of /(x, y, z) by 7), (x, y, z) is: 

Theorem 2.4.3. 

Let f E and y/E S. Then: 


(i) There exists a constant M < max(2,0) , such that 



f < MK2,n^^(f',l/ n) 0,n oo I 


(ii) 

f-Tn 

= 0(n n -> °c, iff / e ; 

(iii) 

f-Tn\ 

1 = ), rt -> oo, iff/is a constant; 

(iv) 

f-T. 

= 0(i/r («-""■ )), n ^ oo, iff K2^^^ (/; r) = 0(i/^(r )), t 0 ; 

(V) 

f-T„ 

= o(ij/(n~^'”)),n-^ oo, iff K 2 ,n;i{f\t) = o(i/r(r^'”)), t -> 0; and, 

(Vi) 

f-T. 

- K2^,^{f•Xln), if S. 


Numerical results corresponding to 2-D Butterworth filter are discussed in 
section 2.5. Segmentations obtained with 2-D Butterworth filter for different values of 
parameters d and m show similar results. Plots of percentage relative errors in the 
approximation with 2-D Butterworth for parameter 0 with respect to different values 
of m show that the smoothing increases with the increase in the value of 6 and 
decreases with the increase in the value of m. 

Chapter III studies segmentations using the discrete convolutions with 
Korovkin, Jackson, generalized Jackson and De La Vallee-Poussin filters. Since their 
application to 2-D and 3-D images is iterative in nature, the degree of approximation 
results considered are 1-D. With the discretizations defined as: 

where (x) = 1/2 + Pk,i(n) cosAx > 0 , x e [0,2;r]. 

For a general following is the asymptotic error formula in the 

approximation of /(x) : 

Theorem 3.2.2. 

If l{n) > m{n) + 2 ,/is a bounded 2;r-periodic function, and f" exists at a point x, 

Kn.m (/’ - /(^) = (1 - Pi.„ )/'(^) + 0(1 - Pi.« )> 

iff for some k>2, 

lim„_« (1 - Pfc.„ )/(l - Pi, J = 



Section 3.3 studies the direct and inverse results on the approximation of 
Korovkin, Jackson and generalized Jackson operators. Let € 2 ^ 1 ^ usual denote the 

space of 2;r-periodic continuous functions on R, and C2^^the sub-space of functions 
having continuous second order derivative. Let / =maxf(x) , /£ and let 

g\2 =81+8", 8^C2^^- Let the Pcctrc s K function3.1 for defined 3.S1 

K{f\t) = M{f -g +r^.g'|:g6 C2;,^},r>0. 

The direct and inverse results on the degree of approximation of f by the 
discretizations are as follows: 

Theorem 3.3.3. 

Let f(x)EC 2 n> and y/eS. If 1-Pi_„ = 0([m(n)]"^) ~ 1/n^ and l(n)>m(n) + 2, 
then: 

(i) - f < AA:(/;(1- Pi where A is a constant; 

(ii) - / = 0(vr(n-2)), n -4 oo, iff K{f-,t) = 0(vr(r")), r 0; 

(hi) =o(y/(n-^)),n^oo, m = 

(iv) ~ K(/;n-'), if S. 

For the De La Vallee-Poussin kernel there holds: 1- Pi,„ = 0(1/ n). Hence the 
above results are not applicable to the discrete filters V,; ;(„)(/; x), which are studied 
separately in section 3.4. 

However, since the relation 1-pj^ = 0([m(n)]~^) ~ 1/n^ is valid for the 
earlier mentioned Korovkin, the generalized Jackson, and Jackson kernels (the case p 
= 2, in the following), for the discrete convolutions A„_2,/(„)(/;x), L„p-p,i{n)U'\^) 

and L2„_2 ,/(„)(/;x), we obtain the following: 

Corollary 3.3.5. 

Let /(x) 6 C 2 n > and yreS. Then, if l{n) > n : 



(0 At- 2 ,i(n)f~f for a constant A; 

(ii) Ai-umI *“ / =0(i//'(n~^)),n->oo, iff a» 2 (/;r) = O(i//-(/^)),r-^0; 

(iii) Ai~-2,l(n)f “ / =o(y/(n ^)),n-^oo, iff a) 2 (/;r) = o(t/A(r^)),r->0; 

(iv) A-uMf-f 

Corollary 3.3,6. 

Let f{x ) e C 2 ^ , p ^ 2, and )//■ e S'. Then, if l{n) >np-p + 2\ 

(i) Kp-p,i{n)f~f ^ A<2)2(/;l/n)-»0,n-^oo, foraconstantA; 

(ii) Kp-p,mf-f =0(y/(n-^)),n-^oo, m(o,if-t) = O(ii/(t^)),t^0- 

(iii) Ap-p,i{n)f-f =o(ii/(n~^)),n-^oo, iff Q)^(f-t) = o(y/(t^)),t-^0 ■ 

(iv) L„p-p,i(n)f-f ~ ^y 2 (/;«‘‘). if ( 02 (f;A^)e S. 

In particular for the Bemsteinian orders n~“, 0< a<2, there hold: 

Corollary 3.3.7. 

Let f{x) e C2„ , 0 < a < 2, and l(n) > n. Then: 

(i) 4-2, /(n)/-/ =0(«-“),n-4oo,iff G)2(/;r) = O(t“),r-40; 

(ii) A„_ 2 J(„)/-/ =o(n"“),n-^oo,iff (y 2 (/;r) = o(r“),r->0; 

(iii) A„_22(„)/~/ ^ n ^ , if ( 02 {f ',t) t . 

Corollary 3.3.8. 

Let f{x) G C 2 j: ,0 < a<2, p>2, and /(n) >np- p + 2. Then: 

(i) Ap-pmf-f =0(n"“),n->oo, iff £U 2(/;0 = O(t“),t-^0; 

(ii) L„p-pMn)f-f =o(n”“),n-^oo, iff W2(/:0 = o(r“),r-^0; 

(iii) ^np-p,l{n)f ~ f if 0)2if\t) ~ t . 



Approximation of De La Vallee-Poussin is studied in Section 3.4, where the 
following results are proved: 

Theorem 3.4.5. (Approximation of derivatives). 

If / is a bounded 2;r-periodic function such that for some natural number m, and at 
some point x, (x) exists, then with l(n) > n+[m/2]+l: 

Moreover, if (x) exists for all x and belongs to C 2 J 1 , the convergence is uniform 
in X. 


Theorem 3.4.6. (Voronoyskaya Type Asymptotic Formula). 

If/ is a bounded 2;r-periodic function such that for a natural number m, at some point 
X, f^"'*'^\x) exists, then with /(n) > n+[m/2]+2; 

Moreover, if (x) exists for all x and belongs to Cin, the convergence is uniform 
in X. 


Let Ca/, (m > 0), denote the space of 2:;r-periodic functions having continuous 
m-th order derivative on the real line. Let |1/|| = max \f{x)\,f e Cjn- Let a Peetre’s K- 
functional, related with the pair {€ 2 ^, 02 ^'^^), for/e 02 /*, be defined as: 

In order to handle the approximation of functions as well as their derivatives 
simultaneously, we would naturally identify /^°^(x)with/(x) , and, C 2 n with C27f 

The following is the main direct and inverse result about approximation of/ by 
Vnn„)f (cast m = 0) and its continuous derivatives /^"“^ by the m-th derivative 

Of Kiwf ■ 

Theorem 3.4.10. 

Let fix) E C2n, m > 0, lin) > n+[m/2]+2, and y/ e S. Then: 



0) f ~ 0, n oo, A being a constant; 

(ii) =0(ii/(n-^)),n^oo, iff = 0(xif(t^)lt -^0; 

(iii) =o(v^(n-‘)),n-^oo, iffK,„^„,^{f-t) = o{w{t^)),t-^fi\ 

(iv) V„,/(„)^"’V-/^'"^ ~ if 

Numerical results on the approximation of Korovkin, Jackson, generalized 
Jackson and De La Vallee-Poussin are presented in sections 3.5 and 3.6. Blurring is 
more on images convolved with Korovkin filter as compared to with that of Jackson 
and Generalized Jackson operators corresponding to the cases n = 128, 64, 32, 16 and 
8. Segmentation result obtained with Korovkin operator for n = 128 seems to be 
practically useful. Spatial convolutions with De La Vallee-Poussin operator for n = 
128, 64, 32, 16 and 8 gives too much blurring, but the result obtained corresponding 
to case n = 4096 is useful for practical purposes. 

In Chapter IV, the segmentation methods based on maximum likelihood 
estimation in univariate and multivariate cases are discussed. The segmentation 
algorithm for the single contrast is given in section 4.2. Section 4.3 contains the 
derivation of iterative scheme and the corresponding algorithm for multivariate 
segmentation (R.K.S. Rathore, S. Datta, R.K. Gupta, S.B. Rao and R. Verma [2001]). 
A detailed study of two pathologies - presence of edema in brain parenchyma, and 
post-traumatic epilepsy with perilesional gliosis, using segmentation methods 
discussed in previous sections, has been done in sections 4.4 and 4.5. Tracing of 
boundaries using cubic splines (S. Datta, R.K. Gupta, S.B. Rao, V.S.N. Kaliprasad 
and R.K.S. Rathore [2000]) is described in section 4.6. 



CHAPTER I : THE GAUSSIAN FILTER AND THE 
MARR-HILDRETH AND RELATED 
APPROXIMATIONS 


This chapter studies the degree of approximation properties of the Marr- 
Hildreth and other related convolutions with the Gaussian function and the resulting 
segmentations. Convolution of an image with Gaussian filter is described in 
introductory section. Section 1.2 contains the basic results on degree of approximation 
by a sequence of linear operators which are required for the subsequent analysis. 
Section 1.3 contains the main result of asymptotic estimates of Marr-Hildreth 
operator. Computation of gradients of arbitrary order has been done in section 1.4 and 
the degree of approximation by Gaussian convolution has been studied in section 1.5. 
Section 1.6 discusses the saturation of convolutions with Gaussian filters and the 
asymptotic error estimates for 2-D Gaussian filters are obtained in section 1.7. 2-D 
gradients of arbitrary order are considered in section l.,8 and the error for smooth 
gradient functions in plane is studied in section 1.9. Finally, section 1.10 includes the 
numerical simulations and discussions for the Gaussian filter and Marr-Hildreth 
operator. 


1.1. INTRODUCTION 

MRI produces 2-D cross-sectional slices containing the 3-D information of the 
examined organs. In the commonly used method where the physician analyzes the 
sequence slice by slice, it is very difficult to reconstruct mentally the 3-D shape of the 
structures explored. To overcome this problem by computer-generated perspective 3- 
D displays of surfaces, the surfaces of the objects to be displayed must be determined. 
Typically, this is done manually slice by slice. Due to their characteristic gray levels, 
bony structures in CT images are amenable to threshold methods. Extraction of the 
surface of soft tissue objects like brain or muscles is much more difficult, as the 



contrast between these objects is poor. In MR images, even though the contrast 
between different tissues is much greater than in CT images, there is not a unique 
correspondence of gray level ranges to different tissues. This necessitates the use of 
an edge detection operator to find the surfaces in the image sequence allowing a 
segmentation of main anatomical constituents an anatomist would like to dissect. As 
discontinuities lead to wrong or misleading structures in the images of the objects, for 
a 3-D reconstruction, it is very important that the surfaces found are closed in 3-D 
space. This requires a segmentation method that recognizes the information in the 
data volume in all three dimensions. 

Following M. Bomans, K-H. Hbhne, U. Tiede, and M. Riemer [1990], among 
various algorithms for edge detection, the Marr-Hildreth operator and related 
approximations: (i) easily extend to 3-D and deliver closed contours without much 
computational effort, (ii) calculate the surfaces based on larger neighborhoods than 
classical edge detection operators, and (iii) work, as claimed in D. Marr and E.C. 
Hildreth [1980], similarly to the low-level processing of the human visual system. The 
2-D Marr-Hildreth operator is defined (E.C. Hildreth [1983]) as: 

C(x,y) = V'(/(x,y)®G(x,y,cr)), 
where I(x, y) is an image, 

G(x, y , a) = (2;7:cr ^ exp[-(xH y ^ ) /(2 c7 ^ )] , 

is the Gaussian function, 

+9^/3y^), 

is the Laplacian, and C(x, y) is the contour image. The convolution (®) in 2-D is 
defined as: 

l{x, y) ® G(x, y, cr) = I (t, s) * G(x -t,y-s,G )dtds . 

• (—00,00) • (—oo^co) 


The Marr-Hildreth operator first smoothens the data volume to remove 
structures of higher frequency components, thereby shortening the frequency range on 
which intensity changes can occur, and only intensity changes at a specific resolution 
are considered. Marr and Hildreth advocate Gaussian operator as an optimal 
compromise between conflicting smoothness constraints as Gaussian distribution is 



smooth and localized both in the spatial and frequency domain. It is thus unlikely that 
artifacts are introduced by the filtering process. Another advantage of the Gaussian 
function is that the amount of smoothing and thereby the considered frequency range 
can easily be adjusted by choosing an appropriate value for the standard deviation. If 
a low standard deviation is taken, only values inside a small neighborhood are 
smoothed; with a greater standard deviation, values over a larger neighborhood are 
smoothed. 

If there is an intensity change with a particular orientation in an image, there 
will usually be a peak in the first, and a zero crossing in the second directional, 
derivative taken perpendicular to the orientation of the intensity change. The detection 
of edges with a directional operator, however, has the disadvantage that the operator 
must be applied at different directions in order to get intensity changes of all 
orientations. Marr and Hildreth chose the Laplacian as a rotationally invariant 
operator, which therefore detects edges of any orientation. The 3-D extension of the 
Marr-Hildreth operator is defined by 

C(x, y, z) = (/ (x, y, z) ® G(x, y, z, a)) = l(x, y, z) ® V^G(x, y, z, cr) , 

G(x, y, z.fT) = (27r(T^)“^'^ exp[-(x^ + y^ + z^)/(2(T^)] , 
={d^/dx^+d^/dy^+d^/dz^), 

V^G(x,y,z,a-) = (2;r)”^^^exp[-(x^ -i-y^ +z^)/(2(7^)][(x^ -t-y^ +z^)/ct^ -3]/cr^ 
As the contour volume can be calculated by convolving the data volume with the 
Laplacian of a Gaussian function, the operator V^G(x, y.z.cr) is called “Laplacian of 
a Gaussian” or V^G -operator. 

The 3-D convolution with V^G is not separable into three one-dimentional 
convolutions. Nevertheless, the 2-D Laplacian is decomposable into two separable 
functions (A. Huertas and G. Medioni [1986]), which for 3-D, results in the 
following; 

V ^G(x, y, z, cr) = G'(x, a)*Giy,<7)* G(z, a) + G(x, a) * G'(y, cr) * G(z, cr) 

+ G(x, cr) * G(y, (7) * G'(z, or) , 


where 



G{^,a) = (2;z:a2)-'^^exp[-(^^)/(2cr2)], 

G"{^,a) = {lna^r^''^{^'^ I a'^ -l)exp[-(^^)/(2a^)]. 

The algorithm needs eight 1-D convolutions and a memory of 3.X.Y.M elements 
(X,Y: size of the data cube; M: filter width) for obtaining the convolution. Also Marr 
and Hildreth have shown that their operator can be approximated by a separable 
operator, the “Difference of Gaussians” (DOG) operator: 

V^Gix,y,z,(J) ~ G(x, y, 1 , 0 - G(x, y, z,a 

» Gix,o^) >!-• G(y,<Te) * Giz,o^) - Gix,o) * GCy.cj;) * G(z,Oi), 
which needs only six 1-D convolutions and approximately 2.X.Y.M memory 
elements. For implementation it has been suggested to prefer the DOG operator for its 
lower memory requirements and smaller number of convolutions (M. Bomans, K-H. 
Hdhne, U. Tiede and M. Riemer [1990]). 

It is reported that if the DOG operator is used as an approximation for the 
V^G operator, Og and O/ have to be determined such that the bandwidth of the filter is 
small and the sensitivity is large: if O/ : Og approximates 1, the DOG and the V^G 
operator are (up to a constant factor) identical, but the sensitivity is 0, and if the ratio 
increases, the sensitivity increases, but the approximation becomes worse and the 
bandwidth of the filter enlarges. As a compromise, Marr and Hildreth chose a value of 
1.6. Some mathematical properties of the Marr-Hildreth operator are studied by V. 
Torre and T. Poggio [1986]. Reportedly the Marr-Hildreth operator detects coarse 
outline of different tissues e.g., skin or ventricular system quite well, but in the region 
of brain sulci the contours found relate to white-gray matter transition than to the 
transition from gray matter to CSF. To overcome the problem it is suggested to shift 
the contours in the sulci outwards using morphological dilation, erosion and closing 
filters for two dimensional objects that alter the form of objects and recognize special 
objects (R.M. Haralick, S.R. Sternberg and X. Zhuang [1987], J. Serra [1982]). 

In view of the fact that it is not very clear as to precisely what smoothing and 
filtering does to the contours, edges and the boundaries associated with different 
segments or objects, we look at the effect of the degree of approximation on the 
segmentations and associated curves. The experiments in subsequent chapters seem to 



confirm that at the present de-facto standard of 256x256 resolution, the segmentation 
properties depend essentially on the degree of approximation of noise reducing basic 
operators which are essentially convolutions based bell shaped functions. 


1 . 2 . BASIC LEMMAS ON DEGREE OF APPROXIMATION 

Some basic results from R.K.S. Rathore [2000], for studying degree of 
approximation by a sequence of linear operators (L„ }, are summarized as follows: 
Let Y cz X he Banach spaces and 0 a semi-norm for which Y = {f e X :<!>(/ )<°°}, 
and let 

Such a A!” is called a Peetre’s /f-functional corresponding to the triple (X, K, O}. 

For a sequence {L„ } of linear operators from X Y, and a sequence } of 

positive numbers such that > ca„ — > 0, n ^ ©o , an inequality of the type 

LJ-f 

is called a Jackson inequality, 

«{LJ) S Sa„-‘ / , 

a Bernstein inequality, and 

^iL,g)<C^(g),gEY, 

the uniform O-boundedness of {L„ } over F. In such a case we say that { F, <I>, AT} is a 
J-B triple with the associated rate sequence {(T„ } . 

Lemma 1.2.1. 

For a continuous non-decreasing function v^(r) satisfying 0-y/(0) <y/{t) for t > 0, 
the following statements are equivalent: 

(i) 

•[ 0 , 1 ] 

(ii) 

• [0,1] 



(iii) For some c> 1 , lim sup^^Q y/(ci)l\j/(t) < c ; 

(iv) For some Ad >0, V^(A5)/i/r(<5) < (1/2)4 (A > A), A<5 < 1/2) ; 

(v) lim^^o supo<,^i /2 v\j/(t)/y/{vt) = 0 ; 

(vi) supo<,<i /2 v\irit)/y/ivt) = 0(1/ 1 Inv |) -^ 0, v -4 0 . 

If y/ is as in Lemma 1.2.1, and satisfies any of the equivalent conditions (i)-(vi) in it, 
we say that y/e S (the class of steady order functions). 

Lemma 1.2.2. 

Let {(T„} satisfy ct„+i > ccr„ — > 0, n oo , K(t)>Q be non-decreasing, and 
(5(r)>0,re (0,1]. Let K((7j)<M{5((jJ + (7jK(Gj/aJ, 1 < n < ; < «, for a 
constant M. If y/e S, then: 

(i) K{t) = 0(y/(t)), if 5(g„) = 0(y/iGj) ; 

(ii) K(t) = o(yr(t)), if <5((T„) = o(i//((jJ); and, 

(iii) K(g„) = 0(5(g„)), if KeS. 

1.3. ASYMPTOTIC ESTIMATES FOR MARR-HILDRETH 
OPERATOR 

Several authors, including K. Yoshida [1971], P.L. Butzer and H. Berens 
[1967], H.S. Shapiro [1969], etc. have studied approximation properties of one 
dimensional convolution with the Gaussian filter: 

G(x,g) = {InG^y^'^ exp[-x ^ /(2cr^)] . 

These studies enable a better understanding of the Marr-Hildreth operator. In what 
follows, we consider the pointwise approximation of gradients of a gerneral order by 
the corresponding gradients of the 3-D Gaussian filter. The main results clarify the 
asymptotic dependence of the degree of approximation on the gradients of higher 
orders as well as the interplay of the smoothness of the gradient in its degree of 



approximation. We begin by developing an asymptotic error estimate and a more 
detailed asymptotic expansion of it for the Marr-Hildreth operator. 

A bounded continuous function f(X) = f{{x,y,z)') is called m-differentiable 

at a point Z e R^, if it admits the following expansion in a neighborhood of the point 
Z: 

/(Z + H) = /(Z) + (H • V)V(X) + o( //'”), H G . 

The function /(Z) is called uniformly m-differentiable on a set 5 c if given an 
arbitrary £ > 0 , there exists a <5 > 0 , independent of Z, such that whenever H <5 , 

f(X+H)-\f(X) + (k\r ^ (H ■ V) V(X)} < £ H . Z 6 S. 

In conformity with = (3^ Idx^ + 3^ / dy^ + 3^ /3z^) , let 

=(3V3x^+3V3y^+3V3z^)^ 

In the sequel n!!, the semi-factorial of n, is the product of all even or odd natural 
numbers, according as n is even or odd. For the 3-D Gaussian filter 

G(Z,<t) = exp[- Z ^ /(2cr^)] , 

there holds the following asymptotic error estimate: 

Theorem 1.3.1. 

If 7(Z) is 4-differentiable at a point Z, 

lim^_oCr-4j'(^)®V"G(Z,(T)-V2/(Z)}=V^/(Z)/2, 
and this relation holds uniformly in Z s 5 if /(Z) is uniformly 4-differentiable in 5. 
If 1(X) is 2p-differentiable at a point Z, 

/(Z)®V^G(Z,<7)-V^/(Z) = 5^^(2A:)a^^*'*^X^_^_^((20!!(2m)!!(2n)!!)'^ 

x(3''=/(Z)/3x2'3y"'"3z''') + o(a'(^-'^), 

where the sums run over 2 < k < p, I, m, n > 0, with l+m+n = k. Moreover, the 
relation holds uniformly in Z e S' if /(Z) is uniformly 2p-differentiable in S. 



Proof: In the following a triple integral [.]dH is an integral over R^, where in each 

» * » 

of the iterated integrals the limits are from - oo to + o® . We have 
/(X)0V2G(Z;a) = /(Z)0((2;rr3'V-^exp{-|x|'/(2(72)}(-3 + |x|'/(72)) 

= ) + Ei.*. 2 p (-^ • V)^ /(Z) + o(\H )] 

x(2;r)-^^V“^exp{- H^/(2 ct^)}(-3+ Hflo^)dH. 
Using (2;r) ^ exp(-r^ /2(T^)dr = 0, ifj is odd and (/-l)!!(y, if j is even 

X AV'”v" exp{-(A^ +/i^ +v^)/(2cr^)}[-3 + (A^ +/i^ +v‘^)l a^]d?ui^v 

X exp{-(A^ +/i^ +v^)/(2cr^)}[-3 + (A^ + +v^)/ G^]d?JiJdv 

= '^''Eos«,E,,„,^„,,„„..[l'((20K2m)!(2n)!)](3=V(X)/32.^3y^"3z^") 

x[(2Z - l)!!(2m - l)!!(2n -l)!!(-3 + (2/ + 1 + 2m + 1 + 2n + i))]( 7 ( 2 '+ 2 m+ 2 n) 

= (3 V3x^ + 3^ /3/ + 3 V3z^);(Z) + X^,,,E,.„,^„j„„..2'^'«2')”(2m)!!(2„)!!) 

X ^ 2(^-1) 0 2<: J )/ ^^2/ 0^ 2m _ 

For the o{H^^) -term, given an arbitrary e > 0, there exists a 5 > 0, such that 
o{H^^) <{e+ I5 ^)H^\He rI Using, 

“ ’ Z G{X , G)dX =0{G^P ^^^ ), (7 ^ 0 , 

the net contribution of the o(H^‘’)-term being o(cr^^^“^^), and the result for 2p- 
differentiable I{X) follows. In particular, for p = 2, 
]im^_^G-^[I{X)®V^G(X-,G)-V^IiX)] 



= 4[{a^ idx^ + aVay" + a" /a^" }/(x)/8 

+ {a'^/ay^az^ +a'^/az^aA:^ +d^ ldx^dy^]l{X)IA] 

= {ll2){d'^ Idx^ +d'^ !dy^ +d^ Idz^f I{X). 

The uniformity of the asymptotic estimates in the cases of uniform differentiability is 
clear as the 5 in above could be chosen to be independent of Z 6 5. 

In the sequel, a similar order of approximation error would be observed for a 
large class of 3-D convolutions with combinations of the Gaussian and other filter 
functions, which at the same time are separable into three 1-D convolutions. In 
particular, for the error in the pointwise approximation by the Gaussian filter; 

Theorem 1.3.2. 

If 1{X) is 2p-differentiable at a point Z, 

2a-2[/(Z)0G(Z,cT)-/(Z)]-V'/(Z) = 2X,cr2^'^-')£^_^_^((2/)!!(2m)!!(2«)!!)-' 

X O^^/CZ) / ) -h o((j^'^-‘^ ) , 

where the sums run over 2 < k < p, I, m, n > 0, with l+m+n = k. This relation holds 
uniformly in Z e 5 if /(Z) is uniformly 2p-differentiable in S. In particular, if /(Z) 
is 4-differentiable at a point Z, 

lim^^oCr"'[2(T-'{/(Z)®G(Z,a)-/(Z)}-V'/(Z)] = V^/(Z)/4, 
and the relation holds uniformly in Z 6 5 if /(Z) is uniformly 4-differentiable in S. 
Proof: Following the proof of the previous theorem: 

/(Z) 0 G(Z,o-) = /(Z) 0 (( 27 r)"^'V"^ expf Z ^ /(2cr^)| 

= '.'t (^) + Ei<.< 2 p • V)^ /(Z) + o(|// )}(2;r)-'' V-' 

xexpf H ^ /(2c7^)}fif. 

As before, the o{H ^^)-term on the right hand side contributes to o(cr^^), which 
holds uniformly in Z 6 5, in the uniform differentiability case. For the remaining term 



.‘■‘teoafl/*!)'' (-« ■ exp{- X ^ K2<r^)]dH 

x_''_'aV”v" exp{-(l^ +v^)K2a^)]dMfidv 
= EoaspS,,,„.,*o.,«,,„.i['/«2')!(2-n)!(2n)!)]0"/(X)/ax“3/"3z^")(2;O-='^CT-= 
x_'_'_’A“ai^"v^" exp{-a’ + 1 ^^ + v’)/(2o-^))<iU/j<iv' 

= 2oa.,L,„«,,„„.»K2()!(2m)!(2n)!)-'0^‘/(X)/3x^'3/”ax^") 
x[(2/ - 1) !!(277 i - 1)!!(2/z - i)!!]o-( 2 /+ 2 «,+ 2 «) 

= Eoax,2..,.«j™«..«2')!!(2m)!!(2„)!!)-'a=‘3«/(X)/3x^'3).^”3z^". 

Hence, 

l(X)®G(X,<7) = IiX) + <THd^/dx^+d^/dy^+d^/dz^)IiX)/2 

+ a\d^ /dx^ +d^ /dy^ +d^ ldz^fliX)l^ 

+ 23x»x,‘^“2,„,,„j*„.„.«2')”(2m)!!(2n)!!)-‘3^‘/(X)/3x“3y^”32^"+o(a^'), 

and the result follows. 

Corollary 1.3.3. 

For the ‘Difference of Gaussians’ (DOG) operator of Marr-Hildreth: 

l(X)®[G(X,a,)-G(X,a,)]-cr,\l-r^W^I(X)/2 

= cr/(l-r^)V^/(X)/8 + o(or/),(T, ->0. 

Note that the approximant in the thorem is nothing but cr“^ -times the error in 
the Gaussian filtering which is separable and can be obtained using three 1-D 
convolutions. Moreover, the asymptotic error is ‘/ 2 -times the corresponding error in 
the approximation of the Laplacian by the Marr-Hildreth operator. This formulation 
for the Laplacian gives good results if liX) is a reasonably noise free float 
reconstruction. However, if I{X) is a noisy or a low order digitized version, a further 
filtering with a G(X,aQ) could correct the situation. Since 

I(X)<S)G(X,(Ti)<S)GiX,<7o) = I(X)®G(X,a,), 



where (K. Yoshida [1971, p. 235]) one is lead to another 

formulation for the Difference of Gaussians (DOG) operator: 

[G(Y,£Tj-G(Y,(j,.)]®/(Y) = G(Z,cr,.)®[G(Z,cro)®/(Y)-/(Z)]. 

1.4. COMPUTING GRADIENTS OF ARBITRARY ORDER 

A class /^.(Y) of functions is called an approximation of /(X). if in some 
sense (e.g., pointwise, uniform, U , etc.) /^(X) /(X),o- -> 0. If for all non- 
negative integers r, s, t > 0, Ig{X)ldx''dy^dz^ I{X)Idx''dy^dz‘ , as 

(T — > 0, /^.(X) is said to possess simultaneous approximation property of an arbitrary 

order. We show in the following that a convolution with the gradients (a linear 
combination of partial derivatives of certain orders) of the Gaussian filter function 
provides pointwise approximation to the corresponding gradients of /(X) of a similar 
order of accuracy as obtained in the previous approximations. The following result is 
crucially needed in proving that (X) = /(X) ® G(X, cr) does possess the pointwise 
simultaneous approximation property of an arbitrary order: 

Lemma 1.4.1. 

For A:> 1, 

a^G(x,cr)/ax* = /(7 ')G(x,o-) , 

where are constants and the sum runs over i, j > 0, such that i-j < k. 

Proof: The result is true for ^ = 1, for, dG(x,G)ldx = -(x/G^)G(x,G) . Hence, 
assuming the result for k, 

a*+^G(x,cr)/ax*^^ = /(j‘^^)Gix,G) , 

which is of the given type. 



Theorem 1.4.2. 

Let /(X)be a bounded continuous function, ^-differentiable at a point X If Q{u,v, w) 
is any polynomial of total degree q in u, v, w, then with Q{D) denoting the 
differential operator Q(d/dx,d/ dy, d/dz), 

lim^_o 1{X) ® [QiD)GiX,<j)] = Q{D)1{X ) , 
the relation holding uniformly in X £ S if I(X) is uniformly <j-differentiable in S. 
Moreover, if 1(X) is 2/7+^-differentiable atX, then 

Qmi(X) ® G{X,a) - /(X)] = Ea.^,v ((2A)!!(2/r)!!(2v)!!)-* 

X (a [Q{D)1{X )J / dx dy-^‘ dz o(a -0 , 

where the sums run over 2< k<p, A,/i,v >0, with 1 + pi + v = k, and the relation 
holds uniformly in X £ 5 if /(X) is uniformly 2p+g-differentiable in S. 

Proof: It is enough to consider the case Q(D) = a'+"'+” ldx‘dy'’'dz" , where l,m,n> 0, 
and l+m+n < q. Since (j-differentiability implies 1+m+n-differentiability of /(X) at X, 

l{X)®Q{D)G{X,(j)='['j{X-H)Q{D)G{H,a)dH 

= ” V(X) + (^!)'' (-H • V)"' /(X) +o( H )] Q{D)G{H,u)dH. 

Using integration by parts, 

• x’^^i^^^^G{X,a)ldx^dy”'dz''dx^- kx'^-^d^^'”^''-^G{X,G)ldx^~%'"dz''dx, 

• /? *^R 

as the integrated part equals zero. It follows that 


0,if ^ < I, 


'd 


l+m+n / 


dx‘dy'^dz" 


a'"+"G(X,(T) 


dy'^dz" 


dx = l\G^’"^iy,a)G’'"\z,(J), 


if jt = l. 


k\ < 


(_ 1 )' x’^-^G{x,a)dx G^”'\y,o)G^'^\z,G), 


if k>l. 


Further, if fc > /, 



•oo 


X 


k 


m—oo 


a'“G(x,cT) ^ 
— dx 

dx%V 


'0,if (^-/) is odd, 


Similar relations are valid with x, >>, z cyclically interchanged. It follows that 
1{X)® Q(D)GiX,CJ) = {il + m + n)\)-H\m\n\[(l + m + n)\/{l\m\n\)]Q{D)I(X) 

+ 3o( )] Q{D)G{H,a)dH 

= Q{D)l{X) + ['\o{H'^"'^’')]QiD)G{H,a)dH. 

Now, given an arbitrary e > 0 , there exists a 5 > 0 , such that 

o( H '*"”)< (e + H " H e R^ 


Hence, taking H = (^,rj,gy, we have from the Lemma that 


'[0(H ‘*"'*")]Q(D)GlH,(y)dH 








. /+m+/z~(r-^)~(r-w)-( v-w) 


•u<m, v,w>0,v-w</i 


) 


_l_ QtV § -2jj. 2+l+m+n-(,r-s)-it-u)-(v-w) 'j 

v^r,^>0,r~j</, t ,u>0,t-u<m, v,w>0,v-w<n / 


<£0(1) + 0(<j^/(5^). 

Due to the arbitrariness of e > 0 , it follows that the contribution of the o-term goes to 
zero as a — > 0, and the required limit relation 

lim^_o /(X)® [2(D)G(X,a)] = Q{D)1(X) 

follows. 

If I(X) is uniformly (j'-differentiable in S, the [£0(1) + 0((7^ term in the above 

is independent of X e S, and the uniformity of the limit relation follows. 

If 1(X) is 2p+(3'-differentiable atX, we have 
Q(D)[IiX)®G(X,(y)-I(X)] 



= e(D)/(X) + ' V)"^ /(X) + o( H )] Q(D)G(H ,a)dH . 


We note that for A: > / + m + n, and r + s + t = k, 

«• % % 

^(_l)r+i«[^!^!^!/{(^_^)!!(^_^)!!(^_^)!!}]^,+,„_(A+^+v)^ 
provided r — X, s — fj,,t — v are all even and non-negative integers, and a vacuous 
semi-factorial is to be interpreted as 1, which arises whenever any of 
r-X,s-/J,,t-v is zero. 

Also, if any of r - A, i' - - v is odd or negative 

'"x''y^z‘[d^''^^''G(X,(T)/dx^dy'^dz'']dX =0. 

% « « 

Hence, 


:;z 


(*!)-' K-H ■V)‘/(X)] Q(D)0(H, a)dH 


JLJl’\-m+n+\<k<q+2p 

= 2, 2 , (*=!)■' E, H(' + 2A)!{m+ 2M)!(n + 2v)!l 
X [(2A) ! ! (2^) ! ! (2v ) ! ! ] “ ‘ }o- 

x[^!/{(/ + 2A)!(m + 2;U)!(n + 2v)!}]32<^^^-"’^^(2(D)/(A:)/a;c^^3y^^3z^'' 

That the contribution of the o( ) -term is o(cr ) is clear from: 

^'”(o(// ^^'”^’'))Q{D)G{H,o)dH 

< ;;;[(£ ^H^I5^)H Ydr,s>.0,r-s<Jr}\^^ lo^) 

< fiOf'V ^q+2p-{r-s)-it-u)-(v-w) j 

^ V^r,j>0,r-j</, t,u>0,t-u^my v,w>0,v-vv<n / 




■u<m, v,w>0,v-w<n 


g-2^2+q+2p-(,r-s)-^t-u)-iv-w) j 


<s0ia^P) + 0(a^’’'^^/8^), 



since 2p+q-{l+m+n)+l-{r-s)+m-{t-u)+n-{v-w)> 2p+q-{l+m+n)>0. If 1{X) is uniformly 
2p+^-differentiable in S, the uniformity assertion for X 6 S follows as in this case 8 
can be chosen to be independent of X 6 S, in which case the [eO((y^^) 
+ / 8^)] - term does not depend on X € S. This completes the proof. 


1.5. ERROR ANALYSIS FOR SMOOTH GRADIENT FUNCTIONS 

Quite an insight into the approximations of gradient functions could be 
obtained by assuming that the images under consideration possess sufficiently smooth 
gradient functions. In fact the images obtained in MRI using FFT on finite raw data 
could be considered infinitely smooth - they being trigonometric polynomials in the 
variables x, y and z- In this section we obtain some inverse results on their degree of 
approximation, which for a bounded continuous gradient Q(D) = Q(d/dx,d/dy, 
d /dz)I(X) to a large extent, characterizes the smoothness of the gradient. Some 
definitions needed follow; Let C denote the Banach space of all bounded continuous 
functions on normed by 

I =sup^^^3{/(X)}, 

and let denote the Banach space of bounded 2-differentiable functions /(X) on 

R^, having all partial derivatives upto order 2 bounded on R^. The space could be 
normed by I = I + I where we define 

|/ 2 = max{||/„|j, .||/^ |,(4/;r) /„ ,(4/;r)||/„|,(4/K) ) 

is a semi-norm on C^. Here the constants (4/7t) have been chosen so as to simplify 
some expressions later. 

Associated with the spaces C and , we define the Peetre’s X-functional by: 
Kil(Xy,t)=^inf^^^,{ l-g +t^ g(X) J/(X)gC. 

It may be noted that for / e C , 

/(X)®G(X,a)-/(X) < "\l{X -aH)-I(X)]G{H,V)dH . 



Hence, if 1{X) is continuous on R^, the integrand being dominated by 
2I G(H,1) and tending to zero pointwise as <t— > 0, it follows by Lebesgue’s 
dominated convergence theorem that I(X)®G(X,a) I(X),(J -^0 . Further, this 
pointwise convergence is easily seen to be uniform on each compact subset of 
(H.S. Shapiro [1969, p. 14]). If, however, I(X) is uniformly continuous on R^. 

co(S) = sup{|/(A: + H) - /(Z)| : |/f I < S}-> 0 , 5 0 , 

and then the integrand being dominated by co(5 H )G(H ,V) , it follows that the 
convergence /(Z) ® G(X,a) I(X), a -> 0, is uniform on R^. 

The degree of approximation of -functions by a Gaussian convolution is 
estimated in the next result, which will also be used in a subsequent result relating the 
degree of approximation with the smoothness of the function being approximated. 

Lemma 1.5.L 

I(X)®GiX,a)-I(X) <3a^ 1\\^,IeC^. 

Proof: If I(X)e , by the mean value theorem of vector calculus, 

I(X +H)- I(X) - H ■ V/(Z) < (l/2)[/zi^/^ -h 

+ 2h2h^ly, +2h^hJ^+2h^h2l^]iX +0H), 

where 0 < 6 < 1 , and all the mixed partial derivatives have been evaluated at the 
indicated point X+0H. It follows that 

I(X)<S>G(X,(J)-I(X) <(l/2)”V^x. 

+ 2^2 I + 2/13 h^\l^ +2hih^ l^]GiH,cr)dH 
<(I 2)' "[hy^ +h2^ +h^'^ +(71 /2)(h2 hs+ h + \ h2)[G(H,a)dH 

<31 ^(t\ 

where we have used the following evaluations (R.K.S. Rathore [1973, p. 65]): 

(271)'^^^ (7~^ t /(2(7^)}dt = (2/7Zy^^(7 , 


Hence, if /(X) is continuous on R^, the integrand being dominated by 
21 G{H,V) and tending to zero pointwise as cr— >0, it follows by Lebesgue’s 
dominated convergence theorem that /(Z)®G(Z,<t) /(Z),cr -^0. Further, this 

pointwise convergence is easily seen to be uniform on each compact subset of 
(H.S. Shapiro [1969, p. 14]). If, however, 1{X) is uniformly continuous on R^. 

£0(<5) = sup{ 1{X +H)- /(X)| : H < S}-^ 0, <5 0, 

and then the integrand being dominated by co(6 H )GiH ,1) , it follows that the 
convergence /(Z) ® G(X,a) I(X), cr 0, is uniform on R^. 

The degree of approximation of -functions by a Gaussian convolution is 
estimated in the next result, which will also be used in a subsequent result relating the 
degree of approximation with the smoothness of the function being approximated. 

Lemma 1.5.1. 

|/(Z)®G(Z,ct)-/(Z) <3(7^ 1\,IeC^. 

Proof: If /(Z) e , by the mean value theorem of vector calculus, 

I(X + H)- /(Z) - H ■ V/(Z) < (l/2)[/zi^/^ + 

+ 2/12^3/^, + 2 h^h^I^+ 2 h^h 2 l^]{X +0H), 

where 0 < 6 < 1, and all the mixed partial derivatives have been evaluated at the 
indicated point Z+ dH. It follows that 

/(Z)®G(Z,a)-/(Z) <(l/2)"V^xc 

+ 2/12 ly^ +2h^ hi +2^ /i2 I ,a)dH 
<U h^+\h hi + h2)\G{H ,a)dH 

<3 / ^ct", 

where we have used the following evaluations (R.K.S. Rathore [1973, p. 65]): 
{2ny^''^(j~^ ' rexpH^ /(2cr^)}dr = {2lnf''^a , 


and 


«* 

This completes the proof of the Lemma. 

Lemma 1.5.2. 

I(X)®G(X,a)\\^<2a-^ l\\,IeC . 

Proof: Using the mles for differentiating a convolution, 

[7(Z)®G(A:,o-)]^ = '''[IiX-H)G^(H,a)]dH 

= '^\l(X-H){-l/a^ +(hy<7^f}GiH,(j)]dH 

<2il\\a-\ 

Also, 

[l(X)®G(X,a)]^ = ‘^[[l{X-H)G^{H,a)dH 

= "\l{X-H)[xyla^}G{H,a)]dH 
<{2l7i)l\a~^. 

Using symmetry, also we have: 

[/(Z)®G(Z,ct)]^^ , [l{X)®G{X,a)]^ <21 a~\ and 

[/(Z)®G(Z,(t)]^^ , [l{X)®G{X,<y)\^ <{2ln)'l a-'^ . 

Hence, 

I{X)®G{X,a) 2 <max{2/ cr~^,{8/7i^) I (j~^]=2a~^ I , 
completing the proof of the Lemma. 

An estimate of error in approximation by a Gaussian convolution in terms of 
the Peetre’s AT-functional is given in the next result: 

Lemma 1.5,3. 

I(X)® G(Z, (T) - I(X) < 3m\a), leC. 

Proof: Let I(X) e C . Then, for g g , we have: 



I(X)®G(X,a)-I(X) < I(X)-g(X) + g(X)-g(X)®G(X,Gy 

+ [g(X)-IiX)]®G(X,C7) 
<2\\l(X)-g(X)\\ + 3a^gl. 

Taking infimum over geC^ we have 1{X)®G{X,g)-I{X) <3K{l\a\ 
completing the proof. 

Theorem 1.5.4. 

Let I{X) & C, and S. Then; 

(i) I{X)®G{X,G)-l{X)\ = 0{\if{a^)),a -^0, m = 

(ii) 1{X)®G{X,g)-I{X) =o(i/r(<T^)),cr -^0, iff = o(v/(r^)),r-40; and, 

(iii) 1{,X)®G{X,g)- 1{X) ~m\G),if 
Proof: That there holds 

I{X)®G{X,g)-I{X) =O(i/r(a^)),cr-^0, if ii:(/;0 = O(vr(r^)),r-^0, and that 

1{X)®G{X,g)-I{X) =o{\if{G^)),G -^0, if K(I-,i) = o(y/(t^)),t-^0, follows 
from the previous Lemma. 

Conversely, if I{X)®G(X,g)- liX) = 0(v/'((T^)),cr — > 0, then 

K{I-t) < K([I(X)- I{X) 0 G(X,G)];t) + mi(X)- g(X)] 0 G(X,(T);r) 

+ ii:(g(Z)0G(X,<7);r) * 

<0(iifiG^)) + 2(tlGf I-g +t^ g^. 

Taking infimum over ge C^, K{l\t)<0(}j/{a^)) + 2{t ! g)^ K{r,G) , and the rest of 
the results follow from Lemma 1.2.2. 

Note that part (iii) of the Theorem justifies the Peetre’s isT-functional error bound 
given in the previous Lemma. The symbol in above, and in the sequel, refers to the 
two quantities being of precisely the same order: there exist positive constants A and 
B such that 

/(X) 0G(X, a) - /(X) <Ais:(/;(T), and, /(X)0 G(X,(t)-/(X) <BX(/;<t), 

for all G sufficiently small. We also note that the familiar Bemstein-Lipschitz orders 
^f{t) = , 0 < a < 2 , belong to the class S. This immediately leads to the following: 



Corollary 1.5.5. 

Let I eC, and 0 < a < 2. Then: 

(i) |/(A:)®G(A:,C7)-/(A:)|| = O(a“),or->0, iff /:(/;/) = 0(/“),r ^ 0; 

(ii) /(X)®G(X,c7)-/(Z)|| = o((t“),£T-» 0, iff ii:(/;r) = o(/“),r-^0; and, 

(iii) \l{X)®G{X,a)-I{X) - m-,a), if if(/;r) ~ 

Let Q{u, V, w) , as before, be a polynomial in u, v, w and let Q{D) = (^dldx, 
didy, d/dz). Then, using l{X)<^[Q{D)G{,X,a)] = QiD)[l{X)®G{X,a)], the 

previous Theorem and the Corollary imply the following characterizations in the 
smooth gradient approximation by convolutions of an image by the corresponding 
gradients of the Gaussian function; 

Jheorem 1.5.6. 

Let Q{D)1{X)e C mdy/eS. Then: 

(i) l(X)®[QiD)GiX,a)]-QiD)I{X) =0iyf(a^)),(7 -^0, 

iff K(QiD)I(X)-,t) = Om^)),.l ^0; 

(ii) l(X)®[Q(D)G(X,a)]-Q(D)I(X) = o(i/r((7")),cr 0, 

iff K{Q{D)I{X)-,t) = 0 ; and, 

(iii) I(X)®[QiD)G(X,G)]-QiD)IiX) ~ m\a), if K{Q{D)l{X)-,t^'^)E S. 

Corollary 1.5.7. 

Let Q(D)I(X) e C, and 0 < a < 2. Then: 

(i) IiX)®[QiD)G(X,cr)]-QiD)I(X)\ = Oia‘=‘),a^O, 

iff K(Q{D)IiXy,t) = O{t‘^),t-^0-, 

(ii) I(X)® [Q(D)G(X,ct)] - Q{D)1{X)\ = o(ct“ ), ct -> 0, 
iff A:(Q(D)/(X);/) = o(r®),r-^ 0;and, 

(iv) l{X)®[QiD)G{X,a)]-Q{P)KX) ~ K{Q{D)l{Xyt) ~ r“. 



1.6. SATURATION OF CONVOLUTIONS WITH GAUSSIAN 
FILTERS 


Finally, we complete the saturation theory for the convolution approximation 
with the Gaussian filter. Let the radial average Ip{X) of a function 1{X)eC be 
defined as: 

/p (X) = {Anp^y^ " I(T)dA{T), p > 0 , 

where A(X,p) denotes the surface of the sphere about the point X of radius p, T is a 
general point and dA(T) stands for the area element on the surface about the point T. 
We call I(X) harmonic if I p(X) = I(X) ,X e R^, p>0. 

Theorem 1.6.1. 

Let /(X)e C. Then; 

(i) I(X)®G(X,<7)-I(X) =O(cr2),a-^0,iff /,(X)-/(X)| =O(r2),r->0; 
and, 

(ii) I(X)®G(X,(y)-l(X) = o((T^),cr^0, iff /(X) is harmonic. 

As an immediate consequence of this Theorem, we have the following: 

Corollary 1.6.2. 

Let Q(D)I(X)g C. Then: 

(i) |/(X) ® [Q(D)G(X,(y)] - Q(D)liX) = 0(ct" ), cr 0, 
iff Q(D)/,(X)-e(D)/(X) =O(r^),r-^0;and, 

(ii) I(X)®[Q(D)G{X,(7)]-QiD)IiX) =o(a^),a-^0, 
iff Q(D)1(X) is harmonic. 

Thus the convolutions with G(X,cr) and its gradients are saturated with the 
order cr^, the saturation classes being functions and gradients uniformly 
approximable by their radial averages with order , and the trivial classes being the 



classes of harmonic functions and harmonic gradients. A proof of the Theorem needs 
the following; 

Lemma 1.6.3. 

If f(X) possesses continuous second order partial derivatives on R^, and, if 
V^/ e C, then: 

fp-f <p^'v^f J6,p>0, 

where equality holds if f{X)= X'^ . 

Proof: We have, 

fp{X)-f{X) = {Anp^y^'\ [f{T)-f{X)}dA{T) 

= (4^)-'* ra/(r)/3n]W)V 

By the Gauss’ divergence theorem, 

/,(X)-/(z)=(4;t)-‘ ' Jr-^ ■■■ rvV(r)MV(r)V 

where S(X,r) denotes the sphere in with center at X and radius r, and dV(T) 
stands for the volume element. Hence 

= p\W^f{X)l6, 
completing the proof of the Lemma. 

Proof of the Theorem: Let us first assume that lyX)- I(X) =0(r^), r— >0. 
Then, as /(A)e C, lyX)-I(X) <Mr^,r>0, for some M > 0. Using Fubini’s 
theorem, 

I(X) ® G(X, cr) -I(X)= ^'\l(X -T)- I{X)]G{T, a)dT 

= ’ U{X-T)-l{X)-\G{T,a)dA{T)^dp 

-[0,oo) • ''A(O.p) J 

= ■ ATip\lJX)-l{X)'\(27ta^)-^‘^&x^{-p^l2a^]dp. 



Hence, 


I(X)®G(X,<y)-l(X) <M 2{p‘^ / G^){2no^y^^^ I2a^]dp 

•[OjOo) 

= 3Mo-^ =0(a^). 

Conversely, if /(Z)®G(X,cr)-/(Z)| =0(o'^), writing: 

/‘"k^) = /(X)®G(Z,T), 

wehave, =0 (ct^). Hence, for some 5 >0, 

l^^\X)®G{X,(y)-l^^\X) <Ba'^, 

for all a sufficiently small. By the asymptotic formula for the Gaussian filter 
/f^’(^)®G(Z,a)-/^^’(Z) = (aV2)V^/t^^(Z) + o(cT^). 

Hence, <2B, and by the previous Lemma, 

vV‘^k^)/6<5pV3. 

Taking limit as T 0, we get: lp{X)-l{X) <Bp^ I3 = 0(p^), which proves part 
(i). For (ii), if I(X) is harmonic, 0 = |./(X)®G(Z,a)-/(X) =o((7^). 

Conversely, from the latter, /^^’(X)®G(Z,cr)-/^’^^(Z) =o((T^), uniformly in t. 
Hence, given any £ > 0, there is a (Tq > 0 such that 

I^^\X)®G(X,g)-I^^\X) <£a^O<(T<(7o. 

It follows that: (X) <2e , and hence by the previous Lemma that 

(/'^’)p(X)-/f'^’(X) <Ep^/3. 

Letting t->0, I p(X) - I(X) <ep^ /3, p > 0. As s is arbitrary Ip{X) = l{X), 
completing the proof 

1 .7. 2-D CONVOLUTIONS WITH THE GAUSSIAN FUNCTION 

Results similar to the 3-D case remain valid for the 2-D case of convolutions 
with the Gaussian function. In this section we consider such results. In view of the 



detailed proofs in the 3-D case, here we mostly outline the proofs to point out the 
differences, if any. In the 2-D case, X = {x, y)' denotes a point in 

^ = [i{d/dx + jid/dy)lV^f(X) = id^f/dx^+d^fJdy^), and 

V^f(X) = (d^f/dx^+d^f/dy^f. 

A bounded continuous function f(X) for which: 

is called m-differentiable at a point X e R^; /(X) is called uniformly m- 
differentiable on a set 5 c R^, if given an arbitrary £ > 0 , there exists a <5 > 0 , 
independent of X, such that whenever H <5, 

/(x+a)->(x)+ (kiy'(.H-v)‘f(X)_ <sH".xsS. 


For the 2-D Gaussian filter 

G(X,<7) = exp[-X ^ /(2cr^)] = (27ra^)"^ exp[-(x^ y^)/(2(y ^)] , 

there holds the following asymptotic error estimate: 

Theorem 1.7.1. 

If /(X) is 4-differentiable at a point X, 

lim^^of^“H^(^)®V^G(X,CT)-V^/(X)}=V^/(X)/2, 
and this relation holds uniformly in X e 5 c R^, if /(X) is uniformly 4-differentiable 
in 5. If /(X) is 2p-differentiable at a point X, 

/(X) ® v 2G(X, tr) - V"/(X) = 5^^ ((2Z)!!(2m)!!)-' 

X 0 / (X ) / dy + o(cr 

where the sums run over 2<k<p, I, m, n>0, with l+m+n = k. Moreover, the relation 
holds uniformly in X e 5 if /(X) is uniformly 2p-differentiable in S. 

Proof: With the double integral [.]dH representing an integral over R^, we have 
/(X)®V^G(X,<t) = /(X)®[(2;r)"^(T''‘exp{-lX ^/(2o-^))(-2+ X ^ /<T^)] 



= + (/:!)-!(-//. V)'^/(Z) + 0( if "")](2;zr)-V-" 

X expj- \h'^ /(2(T ^ )](-2 + \H\/a'^)dH. 

Using expH^ I2a^)dt = 0, if; is odd and (/-l)!!c7', if; is even, 

.*.'(So<*< 2 p • V)* l{X)\2Ti)-^a-^ exp[- \H^ /(2c7^)}-2 + \h\^ lG^)dH 

X AV exp { -(A^ + 2 ) /(2 (t 2 ) } [-2 + (A^ + 2 ) / ct 2 ]^Ai/i 

= Eo...pS/...o./..=.tl/((2/)!(2m)!)]a^*/(X)/ax^'a/'”(2;r)-^ 

x^^ A^V^"” exp{-(A^ +;i2)/(2cr^)}[-2 + (A^ +;i^)/cr^]tfA<i/i 
= Eis«pS,,„K,j™.,2^/«2')!K2m)!!)a»*”-‘>3“7(X)/a^”ay^" 

= (a^/3A:H3V%.^);(X)H-X,a,,S,,,^„_„„.,2t/((2;)!!(2m)!!) 

X(j.2(<:-l)a2*^(-5^)/a^2/a:^2m_ 

For the o(H^^) -term, given an arbitrary e > 0, there exists a <5 > 0, such that 
oiH^^) <{£+ I5^)H^'’ , He rI Using, 

■* Z G{X,o)dX = a -> 0, 


the net contribution of the o( iiT ) -term being o(<7^^^ and the result for 2p- 
differentiable/(Z) follows. In particular, for p = 2, 
lim^_oCT-2[/(Z)®V"G(Z,<T)-V2/(Z)] 

= 4S,,„,o./..=2((2/)!!(2m)!!)-‘a^^/(Z)/ax^'ay^ 

= 4[{avax^+aVay'^}/(z)/8+{aVay^az^+aVa2^ax^}/(z)/4] 

= (l/2)(a^ Idx^+d^ Idy^fliX). 


The uniformity of the asymptotic estimates in the cases of uniform differentiability is 
clear as the <5 in above could be chosen to be ^ dependent of Z € S. 









Theorem 1.7.2. 

If I{X) is 2p-differentiable at a point Z, 

[/(X) ® G{X,a) - /(Z)] = ((T^ /2)v2/(Z) 

the sums being over 2 < k < p, with I, m > 0, and l+m = k. This relation holds 
uniformly in Z e 5 c if/(Z) is uniformly 2p-differentiable in S. In particular, if 
/(Z) is 4-differentiable at a point Z, 

lim^^o<7"H2c^~^[^(^)®V^OT,ct)-/(Z)]-V^/(Z)}=V^/(Z)/4, 

and it holds uniformly in Z £ 5, if /(Z) is uniformly 4-differentiable in S. 

Proof: 

/(Z) (8) G(Z;(t) = /(Z) ® [{2ny^(7'^ exp{- Z ^ /(2a^)}) 

= + )}27r)-' cr-’- 

xexp \-\Hf /(2a^)}iH. 

As before, the o(|jEf ) -term on the right hand side contributes to o(cr^^), which 

holds uniformly in Z e S, in the uniform differentiability case. For the remaining 
term: 

X "aV" exp{-(A^ + ix‘)H2>!^))dXdti 

x[(2!-l)!!(2m-l)!!lCT<^'*^"> 

= Eo.«,S,....,™..«20!!(2m)!!)-a-3“/(X)/3x-3y-. 

So, 



I{X)®G{X,a) = I{X) + a\d^ Idx^ +d^ ldy^)I{X)l2 

+ G\d^ Idx^ +d^ ldy^fl{X)l% 

+ S3a.,'^“L,„>o,,«.t«^)”(2m)!!)-'a“/(X)/3x^'ay^”+o(o-^<'), 

and the result follows. 

Corollary 1.7.3. 

For the difference of Gaussians (DOG) in two variables: 

l{X)®[G{X,G,)-GiX,ai)]-a,\l-r'^)V^l{X)l2 

■ =cr/(l-r^)V^/(X)/8 + o(cr/),cT, ^0. 


1 . 8 . COMPUTING 2-D GRADIENTS OF ARBITRARY ORDER 

Let Q{u,v) be a polynomial of total degree g in m, v and Q{D) = Q{dldx,dldy). 
The 2-D convolution I^{X) = I{X)®G{X,a) also has the pointwise simultaneous 
approximation property of an arbitrary order: 

Theorem 1.8.1. 

If /(Z) is a bounded continuous function, ^-differentiable at a point X, 
lim^_o l{X)®[Q{D)G{X,a)] = Q{D)I{X), 

the relation holding uniformly in Z e 5 c if /(Z) is uniformly (^-differentiable in 
S. Moreover, if /(Z) is 2p-i-^-differentiable atZ, then 

2(D)[/(Z)®G(Z,(T)-/(Z)] = X,<7""j;A.;(((2A)!!(2Ai)!!)-' 

x(a^^[!2(D)/(Z)]/8x^^ay^^)-ho(cr2^), 

where the sums run over 2<k<p, A,/t>0, with X + p, = k , and the relation holds 
uniformly in Z e 5 c R^, if /(Z) is uniformly 2p+g-differentiable in S. 

Proof: It is enough to consider the case QiD) = 9^'*''" /dx^dy'", where I, m>0, and 
l+m < q. Since ^-differentiability implies Z-i-m-differentiability, 


I(X) 0 Q(D)G(X,cr) = ^^I(X- H)Q{D)G{H ,G)dH 

= ' \l{X ) + (/:!)-' {-H ■ Vf I(X) + oiH )] Q{D)G{H,a)dH 

= (/:!)■* /!m![^!/(/!m!)]2(r>)/(Z)+ '\o{H^*"')'\Q{D)G{H ,G)dH 

= Q{D)I{X) + //[o(|h|'^'" )] Q{D)G{H,G)dH. 

As for an arbitrary £ > 0, there exists a 5 > 0, such that 

o(h'^'” <{^+ ,He 


from the earlier Lemma, 


"[0{H'^")mD)G{H,G)dH 


< 


^ - ‘^r,s>0yr-s<i 


X 




.y</ 


'U<m 


/ff') 

‘”’(t)“/Cf')G(H,ff)dH 


< go(y ^l+m-{r-sMt-u) ) 

\^r,s>0,r-s<l,t,u>0,t-u<m / 

+ 0(y §-2^2+l+m-(r-s)-{t-u)\ 

V^r,s>0,r-5</,/,w>0,r-M<m / 


<e0(l) + 0(aV(5^). 

Hence the o-term goes to zero as a -> 0, and /(Z) 0[<2(I>)G(A’,cr)] Q(D)I{X) 
follows. If /(Z) is uniformly ^-differentiable in S, the [£0(1 ) -i-0(ct^ /< 5^)]- term is 
independent of X £ S, and the uniformity of the limit relation follows. 

If /(Z) is 2p+;9-differentiable atZ, 


2(£>)[/(Z)0G(Z,a)-/(Z)] 

= (k!)'‘ (-« ■ V)‘ /(Z) + o( H '*"'■ )_Q(D)G(H,a)dH 

= ■ V)‘ )'Q(D)G(H,c)dH. 

s 



and the o( ) -term contributes to o{(j^^) is clear from: 


/+m , 


[o{H^-"‘m{D)G{H,G)dH 

<"[(e+ H\^ /S^)\H\‘'^^’’]y c f'VfVa") y c ‘"'’( 77 " /cr') 

xGiH,a)dH 

< eOiy ^q+2p-{r-s)-{t-u) | 

<eO{a^P) + Oia^P^^ 15^), 


^ -2^ 2+q+2p-(,r-s)-(,t-u) 


(as 2p+q-Q.+m)+l-{r-s)+m-{t-u)> lp+q-(l+m)>Qi), the asymptotic relation follows. If 
1{X) is uniformly 2p+^-differentiable in S, 8 can be chosen to be independent of X e 

S, so the [£0(cr^^) +0((T^^'^^/<5^)] - term are independent of X € S, completing 
the proof. 


1 .9. ERROR FOR SMOOTH GRADIENT FUNCTIONS IN 
PLANE 

Here, C denotes the Banach space of all bounded continuous functions on 
normed by: I =sup^^^ 2 {/(X)}, and denotes the space of bounded 2- 
differentiable functions I{X) on R^, having all partial derivatives upto order 2 
bounded on R^. The space is normed by ^ C 2 ~ ^ ^ 2 ' where 

/ 2 = max' 7^ , lyy ,(4/71) } is a semi-norm on and 

X(7;0 = inf^gC’^{ 7-g +t^ g(X) 2}7(X)6 C 
defines the Peetre’s X-functional. For 7 e C, 

I{X)®G{X,a)-I{X) = [\l{,X-GH)-I{X)]GiH,l)dH , 
by Lebesgue’s dominated convergence theorem, implies that 



I(X)® G{X,g) I(X), as cr -> 0. 

This convergence is uniform on compact subsets of R^. If 1(X) is uniformly 
continuous on R^, 

(0(5) = sup{ I(X +H)-I(X): H < 5 }^ 0,5 -^0, 
and then the convergence I(X)® G(X,g) -> I(X), cr-^0, is uniform on rI 

Lemma 1.9.1. 

l(X)®G(X,a)-I(X) <(3/2 )ct^ 1 ^,IeC\ 

Proof: If /( X ) e , by mean value theorem 0 

I(X + H)-I(X)-H- V/(X) < (l/2)[/ii^/^ + + 2h,h2l^](X + OH), 

0 < 0< 1, so that 

I(X)®G(X,a)-I(X) <(l/2M/zi^/^ +h2^Iyy +2h^ I^]G(H,(y)dH 
<(/ 2/2)*‘[/zi^+/i2^+(;r/2)(Ai h2)lG(H,(7)dH 
= (3/2)1 y. 

Lemma 1.9.2. 

I(X)®G(X,a) ^<2a~^ I ,IgC. 

Proof: [/(X)(8 )G(Z,ct)L = " 1(X -H)G^(H,G)dH 

= ‘ \l(X - H){-1 /g^ + (hi /a^)^ ]G(H, a)]dH 

m m 

<2I\g~^. 

Similarly, 

[/(X)®‘G(X,a)]3,3, <2Ia-\ 

Also, 

[/(Z)®G(X,£T)]^ = ''l(X-H)G^(H,a)dH 

= '\l(X-H){xyla'^]G(H,a)]dH 
<(2ln)la~'^. 



Hence, /(Z) (8) G(X,g) ^ < max{2 / |a■^(8/:7^^)||/ \a~'^ }= 2g~^ |/||, completing the 
proof of the Lemma. 

Lemma 1.9.3. 

||/(Z) ® G{X,g) - 1(X) < 2K(I-,g\ I e C. 

Proof: Let I(X)e C. Then, for ge C^, we have: 

1{X)®G{X,g)-I{X) < liX)-g(X) + giX)-g(X)<S>G(X,G) 

+ [g(X)-I(X)]®G(X,G) 
<2I(X)-giX) +i3/2)G^ g 

Taking infimum over geC'^, we have I{X)®G{X,g)-1{X) <2K{1-,g), 
completing the proof. 

Theorem 1.9.4. 

Let I(X)e C and y/G S. Then; 

(i) /(Z)®G(Z,o-)-7(Z) =O(i/r((7^)),cr-^0, iff Kil;t) = 0(y/{t^)), r->0; 

(ii) /(Z)®G(Z,(T)-/(Z)| = o(i/r((T^)),(T->0, iff Z(/;0 = o(vr(r^)), t-^O; and, 

(iii) I(X)®G(X,g)-I(X) ~ K{1\g) , if S. 

Proof: 1{X)®G{X,g)-I{X) = O(i/r(o-^)),ifis:(/;0 = 0(vr(r2)), and /(Z)®G(Z 

,CT)-/(Z)| = o(vr(cr^)), if K{l-,t) = o{yr{t^)) , by the previous Lemma. If 

1{X)®G{X,g)- 1{X) =0(y/(G^)), then 

< Ki[I(X) - HX) ® G(X,G)]-,t) + mnX) - g(Z)] ® G(Z,£7);r) 

, +K(giX)®G{X,G)-t) 

<0(il/iG^)) + 2(tlGf l{X)-g{X) +?' g\. 

Hence, taking infimum over geC^ , KiI;t)<0(yf(G^)) + 2(t/G)^K(r,G), and the 
rest of the results follow from Lemma 1.2.2. 



CorollarY 1.9.5. 

Let 1 E C , and 0 < ct < 2. Then: 

(i) I(X)®G(X,a)-I(X) =O(a“),(T-^0, iff /ir(/(Z);r) = O(t“),t-^0; 

(ii) l(X)®G(X,a)-l(X) =o(a“),o--^0, iff Ar(/(Z);0 = oa“),r-> 0 ; and, 

(iii) /(X)®G(Z,cr)-/(X)||~ii:(/;a), if K(l(Xy,t) ~ . 

With Q(u,v) , as before, a polynomial in u, v, and Q{D) = Q{dldx,dldy), the 
previous Theorem and the Corollary imply the following: 

Theorem 1.9.6. 

Let Q{D)1(X)e C, and ij/E S. Then: 

(i) IiX)®[Q(D)G(X,(T)]-Q(D)I(X) =Oiy/(a^)),<7^0, 
iff K(Q(D)l(Xy,t) = 0(y/(t^)),t -^0-, 

(ii) I(X)®[Q(D)G(X,cr)]-Q(D)I(X)=o(y/(a^)),a-^0, 
iff K(QiD)I(Xy,t) = o(y/(t^)),t^O;and, 

(iii) IiX)®[QiD)G{X,G)]-Q(D)I(X) ~K(I;a), if KiQ(D)I(Xy,t^'^)E S. 

Corollary 1.9.7. 

Let Q{D)l{X) £ C, and 0 < a < 2. Then: 

(i) I{.X)®[Q{D)G{X,ay\-Q{D)I{X) =0((t“),ct -^ 0, 
iff ;i:((2(Z))/(.Y);r) = 0(r“ ), t 0 ; 

(ii) l{X)®[Q{D)G{X,<7y\-Q(,D)I{X) =o(ct“),ct-» 0, 
iff K{Q{,D)I{Xyt) = o{t°‘),t-^0\md, 

(iv) I{X)®[Q{D)GiX,oy\-Q{D)I{X) K{r,a), if K{Q{D)l{Xyt) ~ t'=‘ . 



Finally, we complete the 2-D saturation theory for the convolution 
approximation with the Gaussian filter. To this purpose, let the circular average 
Ip (X) of a function I(X) € C be defined as: 

^pU)= \.^I{T)ds{T), 

where C(X,p) is the circle with centre X and radius p, T represents a general point, 
and, ds(T) is the arc length element. We shall call/(Z) harmonic Ip (X) = I(X ) , X 
eR^,P> 0. 


Lemma 1.9.8. 

If I(X) possesses continuous second order partial derivatives on R^, and, if 
V^/(Z)eC, I p(X)- I(X) < p^ 'V^I(x) j4, p > 0, where equality holds if 

liX) ^Xf . 

Proof: 


IpiX)- I{X) = , [liX - pT) - nxws 

= (2;rr'* mX-rT)ldnWdr 

•[0,p]^ *171=1 j 

= {2n)-^' , [V/(r)-n]jAr. 

>|r-X|=r J 


Hence, by the Green’s theorem in the plane, 
IJX)-I{X) = {27t)-^ 


O'* " div{VI{T))dxdy^dr 

.[o.p]^ ..|r-x|<r ^ j 

= ( 27 t)"‘ O'* [V^imidxdy^ dr, 

•'[0,p] ^ * *|r-X|<r j 


so that 


V 2 vT 2 1 




Theorem 1.9.9. 

Let I{X)E C. Then: 

(iii) |/(X)(^G(X,a)-/(X)| =O(a^),a^0,iff|/,(X)-/(Z)| = O(r^),r^0; 


and, 


(iv) I{X)®G(X,a)-I(X) =o(a^), O’ -^0, iff /(Z) is harmonic; 

As an immediate consequence of this Theorem, we have the follow 

Corollary 1.9.10. 

Let Q(D)I(X)EC.Thtn: 

(iii) I(X)®[Q(D)G(X,a)]-QiD)IiX) =0((y^),(7 -^0, 
iff Q(D)I,iX)-Q(D)I(X) =O{r\r-^0-,md, 

(iv) I(X)<S>[Q(D)G(X,<7)]-Q{D)I(X)=oi(T^),(y-^0, 
iff Q(D)I{X) is harmonic. 

Proof of the Theorem: 

If /,(Z)-/(Z) =O(r^),r^0, as /(Z)g C, /,(Z)-/(Z) <Mr^r >0, for an 
M>0. So, 

I(X) ® G(X,cj) - /(Z) = ' '[/(Z -T)- I(X)]G(T, a)dT 

= * [I(X-T)-l{X)]G(T,c7)ds(T)\p 

= ‘ 2np[lp(X)-I(X)](2ncr^)-^exp{-p^f2a^}dp, 

• [0,oo) ^ 

whence 

l{X)®G{X,a)-l{X) ' <{27ta'^y^'^TM ' 2p\2n<7^T^'^ exp{-p^ l2G^}dp 

.[0,oo) 

= 2Ma^ =0(0^). 

Conversely, if /(Z)®G(Z,( 7 )-/(Z) =0(<7^), writing; 

7f"’(Z) = /(Z)®G(Z,T), 

we have, /^^’(Z)® G(Z, 0 ’)-/^^’(Z) = 0(0-^). Hence, for a5 > 0, 

/t^^(Z)®G(Z,<7)-/^^’(Z) <5o■^ 

for all or sufficiently small. By 

/^"^(Z) ® G(Z,(t) - I^^HX) = (cr" /2W^I^'^\X) + oCcr^) , 


we have, (X) <2B , and by the previous Lemma, 

(/'"')p(A:)-/'^'(A:) <p^ /4<5pV2. 

Taking limit as t — > 0, Ip (X) - I(X) < Bp^ 12 = O(p^), which proves part (i). For 
(ii), if I{X) is harmonic, 0= |/(X)®G(X,(T)-/(Z)|| = o(a^). Conversely, from 
the latter, /^^'(X)®G(Z,cr)-/^^'(X) =o((7^), uniformly in t. Hence, given any- s 
> 0, there exists a CTq > 0, such that 

(X) ® G(X,a) - ' < eo•^ 0 < (T < (Tq. 

It follows that: <2e, and 

(I^^^)p(X)-!^^^(X) <sp^l3. 

Letting t^O, I p(X)- I(X) < ep^ /3, p > 0. As £ is arbitrary /p(X) = /(X), 
completing the proof. 


1.10. NUMERICAL SIMULATONS AND DISCUSSION 


In this secton we present some numerical simulation results for the basic 
Gaussian filters with = 128, 64, 32, 16, and 8. 

In Figure 1 . 1 (Original T 2 weighted and corresponding convolved images with 
Gaussian filter), the first row first image is the original cross section, followed by its 
Gaussian filtered version with g’^ = 128. The second and the third row images 
correspond to the Gaussian filter with g' = 64, 32, 16, and 8, respectively in the 
natural order. The filtered image with g * = 128 gives quite a good approximation and 
does not seem to be visibly any different from the original image. The case g'^ = 64 
seems to induce sufficient smoothness without blurring any important structures. 
From g’^ = 32 downwards the increasing blurring is apparent which outlines only 
grosser regions and hides the finer features. 


The percentage relative LI and L2 errors in the approximations for the 
Gaussian filter with g'^ = 128, 64, 32, 16, and 8 are summarized in the Table 1.1, 
which are defined as: 


Percentage Relative LI error = 


^pixels 




^pixels 


/(O 


xlOO, 


Percentage Relative L2 error = xlOO, 


where /(z) and / (z) represent original and filtered images. 


G'^ 

128 

64 

32 

16 

8 

% Relative LI Error 

0.19906 

4.69717 

^9.52364 

14.63921 

' 17.40059 

% Relative L2 Error 

0.17313 

4.23735 

9.62608 

15.87283 

' 19.08742 


Table 1.1: % Relative Errors in the Approximation of the Original Cross-Section with Gaussian Filter 
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Figure 1.2 (MLE Segmented images of Convolutions with Gaussian filter) 
contains MLE-based segmentations associated with the images in Figure 1.1. The 
number of iterations taken for original image is 59400 and for a* = 128, 64, 32, 16, 8 
are 18600, 13550, 3050, 7750, 3700 respectively. The segmentations associated with 
the original and a'* = 128 are unnecessarily detailed. The cases corresponding to = 
64, and cr'^ = 32 seem to be of practical use, while those with a’* = 16, and = 8 
provide grosser simplifications. 

For a"' = 32 and 16, there is a splitting of edema into three regions, whereas 
these merge to give two regions for = 8 and boundary of edema thickens with the 
decreasein the value of o'*. All separate regions corresponding to edema are 
represented as different stages of its growth inside the brain parenchyma. As filtering 
increases, the segments within cranial part merge to give one segment and there is a 
formation of boundary between the major segments comprising of brain and its 
cranial part. Excessive filtering is useful in the separation of brain from its cranial 
part. 


In conclusion, for the Gaussian filter the choice o'* = 64, 32 seem to be 
satisfactory for the purpose of segmentation, and the MLE based segmentation 
appears to be promising. 

Figure 1.3 (Gaussian Difference Magnitudes) presents the original cross 
section of Figure 1.1 and the magnitudes of the differences of the each of the rest of 
the images in Figure 1.1 with the previous one. As these approximate the Laplacian, it 
is clear that the original image possesses non-uniform graininess and that the 
segments of interest do not have uniform intensity along and across the boundaries. 
As expected, the Laplacian boundaries thicken with decreasing a'*. It may be 
concluded that the Laplacian in the present case does not seem to provide closed and 
continuous boundaries. Since filters for the computation of the Laplacian, and the 
DOG filter behave similarly to the absolute differences, their utility for the MRI 
images of the present category may not be very competitive with the MLE based 
segmentation. 



Figure 1.4 (Sobel maps (Gaussian filter)) displays the Sobel maps of the 
images in Figure 1.1. They possess the same general behaviour with respect to 0 as 
the difference magnitudes. However, they provide a much better description of the 
edges as compared with the difference magnitudes. The discontinuities and hon- 
closedness are observable. 

Figure 1.5 (Mathematical Edges of MLE-segmented original and filtered 
images) displays the precise edges from the respective MLE-segmented images. 
Figure 1.6 (Edges traced from original and Gaussian filtered images) displays the 
results of using a certain popular edge tracing commercial program on the filtered 
images. Figure 1.7 (Edges traced from MLE-segmented original and filtered images) 
does the same on the segmented images. 

Results of Marr-Hildreth operator with different values of a (1.0, 1.61, 2.0, 
2.576, 3.0, 4.0) are shown in Figure 1.8. It is observed that the obtained closed 
contours are well represented in the cases of 0 = 2.0, 2.576 and 3.0. Figure 1.9 shows 
the results of Difference-of-Gaussian in two-dimension (A. Huertas and G. Medioni 
[1986]) with the ratio of 1.6 between the values of 0 e = 0.8050 and 0 i = I.2880, the 
different values of 0 are 1.0, 1.2, 1.4, 1.6, 1.8, 2.0. The number of closed contours 
decreases with the increase in the value of 0 , i.e., due to increase in the smoothness 
within an image, which is expected. 
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Figure 1.1: Original T 2 weighted and Corresponding Filtered Images with Gaussian Filter 
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Figure 1.2: MLE Segmented Images of Convolutions with Gaussian Filter 










Figure 1.3: Gaussian Difference Magnitudes 








Figure 1,4: Sobel Maps (Gaussian Filter) 
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Figure 1.5: Mathematical Edges of MLE Segmented Original and Filtered Images 




















Figure 1.8: Edges via Marr-Hildreth Operator for Different Values of a 
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Figure 1.9: Edges via Difference-of-Gaussian Operator for Different Values of g = 0805a,ai = 1.288a) 








CHAPTER II : DIFFERENTIAL SMOOTHING AND 
BUTTERWORTH FILTERS 


In this chapter we study Butterworth and some related filters and specially the 
degree of approximation of function / (x) by trigonometric polynomial filter 

T„(x) = T„(f;x) = <2o(/)/2 + ^j^^^^(a^(/)cosb: + fc^(/)sin bc)/(l + 0(A:/n)^'"), 
where a^{f) and bj^{f) are the Fourier coefficients off, which is characterized by a 

certain minimization property. Section 1.1 is an introduction to simultaneous 
approximations with a minimization of derivatives of different order. The 
approximation properties have been derived in section 2.2. Section 2.3 is a study of 1- 
D Butterworth filters along with its generalization to 2-D, and the 3-D Butterworth 
filters have been studied in section 2.4' Numerical simulations corresponding to the 
approximations and the resulting segmentations by 2-D Butterworth filters have been 
discussed in section 2.5. 

2.1. APPROXIMATIONS ALONG WITH DIFFERENT ORDER 
DERIVATIVE MINIMIZATION 

The behavior of the derivatives of approximations to functions representing 
images could be controlled to quite an extent. The idea is to retain a good degree of 
approximation while controlling the local variations through minimizing the 
magnitude of different orders of derivatives of the functions. This is expected to yield 
better images for the purposes of segmentation as it is well known that the 
segmentation methods do not perform nicely for cross-sections with a high amount of 
local intensity changes. Since the intensity variation is governed by derivatives of 
different orders, it seems natural to try smoothing, without much compromise in the 
approximation of the cross section, while simultaneously minimizing the derivatives. 
Let denote a complex trigonometric polynomial of order n. Since the Bernstein 



inequality implies that the sequence of functions is bounded both in 

2 

^ 2 n as the MRI reconstructions are trigonometric 

polynomials, at least in the x and y-directions, one could derive meaningful separable 
filters by minimizing 


where f{x) represents a complex cross section along rows or columns. HereL 2 ;j^ is 
the space of 2;r-periodic square integrable functions 

/(x) = ag / 2 + (< 2 ^ cos kx->rbi^ sin kx) , 

with the norm given by 


r ' 2 

/ = . ,/W dx 




and 0 is a parameter. If 6 is large the smoothing is expected to get increased. 
Moreover, by increasing m the smoothing effect ought to decrease, as higher 
derivatives have only secondary or finer effects on intensity change in digital images 
(D. Luo [1998]). 


2.2. THE MINIMIZING DIFFERENTIAL FILTER 

Let Tn denote the space of trigonometric polynomials of order n, n> 0. The 
characterization of the form of the trigonometric polynomial Tn minimizing the 
quantity Qm of the previous section is given in the following result. The form of the 
minimizing polynomial makes it depend linearly on the function /. Subsequently the 
approximation properties of r„ have been derived. 

Theorem 2.2.1. 

For f E ,m = 1,2,3,..., and 6 > 0, the trigonometric polynomial Tn(x) e Tn 
minimizing 
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2" /W-r.u) = +0 


is unique and is given by 

W = (/; ^) = <^0 (/) / 2 + (fljt (/) cos b: + (/) sin b:)/ (l + ), 

where (/), bj^ (/) are the Fourier coefficients of the function f{x ) . The map/-^ 
Tn is linear, and, r„ - / 1 -» 0, n -» oo. 

Proof; To find the minimizing , writing 


and using the orthogonality of the functions e‘^. 


^ 2 

The term c^t does not affect the minimization, and it is clear that each of the 

summands (cjt - g^. ^ +6(k/n)^'" must separately get minimized. If = 0, 
gjt =0 is the minimizer. Hence with =^c^., one has to minimize 

g((l>)=l-<l)^+e(k/nf’"(t>'\ 


for which the minimizing 0 is real and is given by: 

0 = g'(<l)) = 2(p{ + eik/nf’^}-l), 

so that 

and as 

g\(l)) = 2(l + e(k/nf'")>0, 

the minimizer is unique. 

Passing to the trigonometric polynomial form, it follows that the. trigonometric 
polynomial r„ (x) minimizing Q„ is given as: 

T„ (x) = T„(f;x) = ao (/) / 2 + (a^ (/) cos kx + b^ (/) sin kx)/ (l + 0 (^ / nf"‘ ), 


with 



= ^ ^ . f(x)coskxdx, bt(f)=n: ^ f(x)sinkxdx, 

•1-71, }r] ' -[-s^] 

being the Fourier coefficients of fix). Finally, for the convergence, 

<'23l\n 

Note that the linear space invariant map T„ : f if',x) can be evaluated by FFT 
and is applicable separably in any of the jr, y and z-directions. Also note that T„ =1. 

We study the effect of this filtering for different 6 and m and also the degree of 
approximation of/(x)by the polynomials T„(f]x ) . Let L 2 „ 2 m, (m > 1) denote the 

subspace of consisting of functions fix) for which 

/2 + 2,,, ( a/ + i., ^ )f ~ . 

With this norm L 2 ;r, 2 m^ ^ Banach space, and we define the associated Peetre’s K- 

functional by 

2 2 

Since all trigonometric polynomials belong to is dense in L 2 „ , from 

which it follows that for all fE.L^J', K 2 „if',t) -^0 , as r— >0. 

Lemma 2.2.2. 

f-T.(f) <maxa.0)n-^"/L. 

Proof: Since, 

fix)-T„if-,x) = ^^^^^^ia,if)coskx + b,if)smkxieik/n)^'"y^^ 

+ ^k>n (/) cos kx + b^ (/) sin kx), 



The main direct and inverse results on the degree of approximation of f{x) by 
T„ (x) are summarized in the following; 

Theorem 2.2.3. 

Let / G L 2 „^, and t// g 5. Then: 

(i) There exists a constant M < max(2, 0) , such that 

/-r„ 

(ii) \f-Tn = ), n oo, iff / g L 2 ^, 2 m^ i 

(iii) ||/ -T„\ = o(n”^'" ), n — > oo, iff /is a constant; 

(iv) f-T„ =0(vA(n-2'")),n-^~, iff (/;/) = Od/r^^-)), t 0 ; 

(v) f-Tn= o(vr(n“^'”)), n -> oo, iff K 2 ^(f;t) = o(y/(t^'")), t 0; and, 

(Vi) f-T„ ~ K^JWn), if S. 

Proof: Since, 

/ (x) - r„ (x) = i^k cos kx + bi, sin kx'jp /[n^'" +6 k^""]) 

+ 'Zk>n cos b: + ^7^ sin kx), 

< /-g +|^-r„(g)| + ||r„(g-/)|<2||/-g +max(l,0)n-2'"|g||2„. 
Taking infimum over g g , we have 

where Mg = max(2,0). As, — > 0, t — > 0, (i) follows. 

If feL 2 „ 2 m, by the Lemma, f -T„\\ = 0(n~^'"),n-^o=. Conversely, assuming 
/ - r„ = 0(n"^"' ) , if Z is a fixed natural number, for all n > / sufficiently large, there 
exists a constant M such that 

^>^2- /-t; 

- Xi<i:<z (f)coskx + b^.if) siu kx)d(k / n)^”’ _1 + 0(^ / 

> '^^^^^j{ak(f)coskx + bi,if)smkxX0k^'")l + 9(k/nf'^ . 



Taking limit as n — > <», we have 

(/) cos kx + bi,(f) sin kx)(e k^”') < M. 

Letting / — > <» , we conclude / e • This completes the proof of (ii). 

For (iii), if / is a constant, f -T„ = 0 = o(n"^'" ) . Conversely, if / - r„ J = ) , 

as in the proof of (ii), we get 

2it>i (/) cos kx + b,^if) sin kx)(e k^’") =0. 

Hence, <3^ (/) = b/^ (/) = 0, for all k >1, and f{x) must be a constant. 

To prove (iv)-(vi), if = 0(r(r^'")), t 0, f -T„ = 0()/r(n"^'")), n -4 oo, 

and if K 2 n(f',t) = o(}j/(t ^’")), t -> 0, |/ -T),! = o(y/(n ~^'”)), follows from (i). 
For the rest, assuming / - T„ = 0(y/(n ~^'^ )), n -> oo, 

K2,„ (/;/)< K2,„ (/ - r„ (/);0 + K2,„ (T„ if - gft) + K2,„ iT„ igft) 

+ f-g +||g 

Taking infimum over g e L 2 n, 2 m > 

(/;r) < 0(v/(n-"'”)) + (m)""- (/;«-')), 
and the results follow from Lemma 1.2.2. 

Taking ^(t) = , for the Bemstein-Lipschitz orders , we have: 

Corollary 2.2.4. 

Let / e L 2 ^ , and 0<a<2m. Then: 

(i) f-T„ =Oin-^),n-^oo,iff K2jf-,t) = Oit^),t^O-, 

(ii) f-T„= o(«'“ ),n-^oo, iff = o(t“ ), r -4 0 ; and, 

(hi) f-r„ ~n-“, if /i: 2 „(/; 0 ~r“- 


Note that in the MRI reconstructions using FFT, the relevant choice 
CT„ = 2"" satisfies (7„^i >(l/2)(7„, so that in the converse parts a given degree of 



approximation along the sequence {2"} itself implies the corresponding smoothness of 
/as given in the theorem and the corollary. 


2.3. THE BUTTERWORTH FILTER AND GENERALIZATIONS 

A lowpass filter produces no or little modification in the Fourier coefficients 
corresponding to frequencies below a given threshold called the cutoff frequency. It 
allows the low frequency components pass and reduces or blocks high frequency 
components. The equations for the transfer functions of the Butterworth low pass 
filters of and the 2"“^ order are given by r(/) = l/(l + (r^/c/^)) and 

T(/) = l/(l + (r'^/c/'^)), respectively, where cf is the cutoff frequency and r is the 
square root of the sum of the squares of the frequency components. The cutoff cf is 
normally in the interval [0, 0.5]. The second order filter is steeper than the first order 
filter. 


The minimizing differential filter of Section 2.2 may be termed one- 
dimensional Butterworth filter of order m. A differential characterization of the two- 

dimensional Butterworth filters may be made as follows. Let Tn ,2 denote the space of 


trigonometric polynomials of order n, n > 0, in the two variables x, y. Let k stand for a 
double index {k^,ky)and \ k +ky^. For members of T„, 2 , then, | k p< n^. We 

will consider the space L 2 „ 2 ^ 2;r-periodic functions square integrable on the 

square {-n,Tt\x[-7t,7t\ with the norm given by 


/ =^. .. .f(x,y)^dxdy^. 


1/2 


Theorem 2.3.1. 

For / e L2^ 2^ = L 2, 3, ..., and 0 > 0, the trigonometric polynomial T„(x,y) e Tn ,2 

minimizing 



01 


”Cj n-»"T,(x.y)ldx-‘dyl "^dxdy. 

is unique and is given by: 

Tn(x,y) = T„(f;x,y) = ^ txp{i(k^x + ky y) }, 

where c^(/) are the complex Fourier coefficients of the function y(x,y). The map:/ 
— > Tn is linear, and, T^- f -> 0, n — > <». 

Proof: Writing 

/Uy) = 2i>o^* exp{z-(^^x + ^yy)}, 
and 

Tn (x, y) = X|i|<„ 8k exp{/(A:^x + ky y) }, 
by the orthogonality of the functions exp{z'(^^x + y) } over [-;r,;r]x[-;r,;r], 

fi- = i.r)+E„>,h ' 

= S|»|S »(^» -Sk ’ +(e/n"")|<:p“ *»/+E|t|>,h "■ 

As ^ does not affect the minimization, each 

must separately be a minimum. If Ck = Q,gk = 0. Else, putting ^ , 

g{^)=\-<l>^+6{\k\lnf^(^^ 

is minimized by a real ^ satisfying 

g'(^i) = 2(4 + 0(|^ I /«)"'” }-l)=0, 

and is uniquely given by: 

0 = _l + 0(|/t|/n)2'"_"\ 

Since = 2(l + 0(| ^ | / n)^'” )> 0 , it is indeed a minima. 

It follows that the trigonometric polynomial Tn{x,y) minimizing Qm is given as: 
T^(x,y)=T„(f-,x,y) = exp{z(^,x + ^ 3 ,y)}/[l + 6>(| k \ /«)''”], 


with 



being the Fourier coefficients of f{x, y). For the convergence, 

ii/-r,ii = ( 2 ^)=is,^, 8 (i i I /n)^" 1 + 9(1 k I +Sw>.hr] 

M2n?W" 1^ + 1^*1' K o. " - 

The linear space invariant map T„ : f ^T^{f-,x, y) has ||r„ | = 1 , and it can be 
evaluated by FFT using multiplicative replacements: -> /[1 + 0 (| yt | /n)^'"] in 

the frequency domain. The effect of this filtering for different 6 and m and also the 
degree of approximation of/(x,y)by the polynomials T^(f-,x,y) is considered 
below. 

2 2 
Let L2„^2,2m . (m ^ 1 ) denote the subspace of L2ir,2 of functions f(x, y) for 

which 

i/L=( 4 - 4 Kr- 2 ;,jkrhrr<»- 

With this norm ||.||2in. ^^271,2,2)^ is a Banach space, and the associated Peetre’s K- 
functional for the pair (L2;j_2^>i^rt,2,2m^) of spaces is given by: 

■^2m,2(/iO =i^f{ f ~ S +?^”' 8 2m' ^n,2,2m \ f ^ ^JC,2 ■ 

Since all trigonometric polynomials in two variables belong to L2„x2m^ > it is dense in 
L2^ 2^ > if follows that for all / e L2„^2^ , ^^2^,2 (/iO as t -> 0 . 

Lemma 2.3.2. 

If /e Vu™'. f-W) S maxa0)n-=^ / „■ 

Proof: We have, 

fix, y) - T„ if; X, y) = {c^ if) exp{i(A:^x + kyy)}) 

x(e(| k 1 /n)^'”)/(l + 0 (r^ i/n)^'”)+ 2 *>„(ci(/)exp{i(M + ^y 3 ')})- 



Hence, 


/ - r, (/) < (/) x»., C/) ^ T 

<max(l,0)n-2"’||/ . 

The main direct and inverse results on the degree of approximation of f(x, y) 
by T„(f;x,y) are: 

Theorem 2.3.3. 

Let / e L2„^2^ , and y/e S. Then: 

(i) There exists a constant M < max(2,0), such that 

f-T„ <M^2m.2(/;l/«)^0,n-^oc; 

(ii) /-r„ =0(n-2'”),n-^oo, m fGL 2 „^ 2 . 2 m^-, 

(iii) f -Tn = oirT ^"' ), n -> oo, iff /is a constant; 

(iv) f-T, =0(v^(n-2'”)),n-^oo, iff is:2,_2(/;0=O(vr(r2-)),t-4 0; 

(V) 7-r„ =o(i/r(n-2'”)),n-^=c, iff ii:2„,2(/;r) = o(v^(r2-)), t 0 ; and, 

(Vi) f-T, ~ /i:2„.2(/;l/n), if 

Proof: Since, 

/(x, y) - T„ (/; X, y) = (c^ (/) exp{/(^^x + ^3, y) }) 

X (0(7 I /n)^"" )/(l + 0(7 I /n)^'” )+ (ci {f)tx^{i{k^x + kyy ) }), 

/-r„(/) < /-^|+ g-r„(g) + r„(g-/) 

<2||/-g +max(l,0)n-^"’ g||2^. 

2 

Taking infimum over g e L2„^2.2m > 

<M,K2„^2(Al/n), 

where Mg = max(2,0) . As, ^2m,2(/ ;0 — ^ 0. t 0 1 (i) follows. If / e L2j^ 2.2m > t>y 
the previous Lemma, | / - r„ || = 0(n“^'” ), n — > <» . 



ij-r 


Conversely, assuming / - r„ = OCn-^"” ) , if / is a fixed natural number, for all n > Z 
sufficiently large, there exists a constant M such that 
M >n^"' f-T„ 

- + 0(1 k I ln)^’\l + d(\ k | 

^ 'Zi<k<ii^k(f)^^V{Kk^x + kyy)}ye\k\'^'")l + e(k/nf"'~^ . 

Taking limit as n — > oo , we have 

Yi<kA^k(f)^^P{Kk^x + kyy)})iek^'”) <M. 

Letting Z -> oo, we conclude that / 6 L 2 „^ 2 , 2 m^- This completes the proof of (ii). 

For (iii), if/ is a constant, / -r„ = 0 = Conversely, if f-T„ = o(n"^'"), 

as in the proof of (ii), we get 

2*>i (/ ) exp {Z(/:_,x + kyy)}\ek ) = 0 . 

Hence, (/) = 0, for all k >1, and / must be a constant. 

To prove (iv)-(vi), if K 2 ^ 2 ^f\t) = 0(vr(r^'")), t 0, \f -T^\ = 0(v^(n‘^'”)), n -> oo, 
and if Ar 2 m, 2(/!0 = o(V'(^^’”)))t-^0, /-/, =o(v/'(:pi"^'”)), n— )-oo, follows from 

(i). For the rest, assuming f -T„ = 0(\i/(n ~^"' )), n -4 «■, 

if 2„,2 if It) < if 2., 2 (/ - T„ (fy,t) + ^2„.2 (T„ if - gy,t) + if2..2 (T„ (^);0 

<0(vr(n-2'”)) + f2-n2'" /_^ + ^1^^ . 

Taking infimum over g e. L 2 „ 2 , 2 m > 

and the results follow from Lemma 1.2.2. 

Taking i/r(r) = ^ Bemstein-Lipschitz orders n"“ , we have; 

Corollary 2.3.4. 

Let / e L 2„^2 > and 0 <a<2m. Then; 

(i) / - T„ = 0(n““ ), n -> o°, iff if 2^.2 if'J) = Oit^),t-^0-, 


(ii) /-r„ =o(«-“),n^oc, iff is: 2 „. 2(/;0 = o(r«),r->0;and, 

(iii) /-r„ ~n-\ if 


In the MRI reconstructions using FhT, the relevant choice ct„ = 2"” satisfies 

^n+i - > SO converse parts a given degree of approximation along 

the sequence {2”} itself implies the corresponding smoothness of/ as given in the 
theorem and the corollary. 


2.4. THE 3-D BUTTERWORTH FILTER 


Following the differential characterization of the 2-D Butterworth filters, a 
differential characterization of 3-D Butterworth filters is considered in this section. 

Let Tn ,3 denote the space of trigonometric polynomials of order n,n>0, in the three 

variables x, y, z. Let k stand for a triple index (k^,ky,k^) and [ k p= kj' +ky^ +k^. 

For members of Tn, 3 , then, \kf‘<n^. We consider the space L 2 «, 3 ^ of 2;r-periodic 

functions square integrable on the cube {-7i,7i\x{-7i,7i\x[-n,n] normedby 


/ = 


r 




* [ -;r ,;r ] • [ -;r ,;r 1 • [ -;r ,7r ] 


2 ^ 
f(x,y,z) dxdydz i 

J 


1/2 


Theorem 2.4.1. 

For / e L 2 ;j 3 ^ , m = 1, 2, 3, . . ., and 0 > 0, the trigonometric polynomial r„ (x, y, z) g 
T n ,3 minimizing 

Q^={2ny^' ' * \f{x,y,z)-T„{x,y,z)^ 

^So</ j (x, y, z)/dx‘dy ^dz^ j^Zx^fydz, 

is unique and is given by: 

T„(x,ya) = T„if;x,y,z) = J,^,^,„c,if)l + d(\k\/nf'"\xmk,x + kyy + k^z)}, 



where e^(/) are the complex Fourier coefficients of the function fix,y,z)- The 

map: /-> T„ is linear, and, T^-f -^0, /2 > CXD. 

Proof: Writing 

/ U 3^, z) = exp{f(^^x + kyy + k^z)}, 

and 

Tn(x,y,z) = &^p{i(k^x + kyy + k,z)}. 

by the orthogonality of the functions exp{i(/:^x + A^y + yfc^z)} over [-: 7 r,;r]x[-;r,;r]x[- 
7r,7r], 

X|A|>n 

= Sws, ( c. - ^ + (e / ) I i P* «/ )+ c. = . 

To minimize Qm, each [^Cp.- gj.^ +0(A:/n)^"' must be a minimum. As gp = 0, if 

Ck = 0, with gp=(l)Cp, 

gi(l>)=l-<t>^+e(\k\/nf’^<l>^ 


is minimized by the real 0 satisfying 

g'(0) = 2(4 + 0(1* I /«)'"’ }-l)=0, 

and is uniquely given by: 

<l)=l + d(\k\fnf’"~\ 

Since g*(<j)) = 2(l -i- 0(| A: | / )> 0, it is indeed a minima. 

Hence the trigonometric polynomial r„ (x, y, z) minimizing <2m is given as: 

T„ (x, y, z) = r„ (f;x, y, z) = ^k exp{i(*^x + y + } /[I + 0(| * | / n)^"* ] , 


with 


c,(/) = (2;r)- 


• [-;r,7r] • [-;r,;r] • [-n;7t] 


f(x, y, z) exp{i(kj^x + kyy + k^z) ]dxdydz , 


being the Fourier coefficients of fix, y, z) ■ For the convergence. 
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f-T,= k \ln)^ 1+9(1 Ic I ,n)\-' g, ^ 

= (, 7 xf^n ”S|j|s^«/+ 2 ni>v;‘^i ^}-» 0 , «-+■». 

As in the 2 -D case, the linear space invariant map T„(f-,x,y,z) has 

= 1 , and can be evaluated by FFT using multiplicative replacements: 

^ ) ill the frequency domain. The approximation properties of 

the polynomials r„ (f;x, y, z) are considered below. 

Let L2nx2m ) (w > 1 ) denote the subspace of ^ of functions f(x,y,z) 
for which 

/ 2m = "F" < “• 

With this norm |l.||2m, l-2n,^,2m is a Banach space, and the associated Peetre’s K- 
functional for the pair (^2^,3^ , spaces is given by: 

2 

Since all trigonometric polynomials in two variables belong to L2„^i^2m > it is dense in 
L2 ;j 3^ , and it follows that for all / e , Ar2m,3 (/; 0 ^ 0 » t -> 0. 

Lemma 2.4.2. 

If /£ i2.,3,2«'. /-r.(/)|<max(l.e)n-^" /L- 

Proof: We have, 

f(x,y,z)-T„if-,x,y,z) 

= Ei<|*|<« (/)exp{i(*^x + + ^^z)})(< 9 (| k \ /«)^'")/(l + 0 (| k | /nf"’ ) 

+ E*>„ (/) exp{z'(^;,x + kyy + k^z)}). 


Hence, 


oo 


/ - r„ (/) < X,.. (<;*(/)')}"' 

<max(l,0)n‘^"'|/ 2 ^. 

The main direct and inverse results on the degree of approximation of 
fix,y,z) by T,^(x,y,z) are: 

Theorem 2.4.3. 

Let fe L 2 ;j 3 ^, and \j/e S. Then: 

(i) There exists a constant M < max(2,0) , such that 

f ~T„ < Mi^’2^_3(/;l/n) — > 0, n — » oo; 

(ii) / - r„ = 0(n-2" ), n oo, iff f e ; 

(iii) f -Tn =o(n~^'”), n->oo, iff/ is a constant; 

(iv) /-r„ =0(i/r(n-2-)),n-^oo, iff is: 2 „_ 3 (/;r) = O(i/r(r2-)),t^0 

(v) /-r„ =o(v/-(n'^'”)),n->oo, iff i5:2^_3(/;r)=o(vr(r^"')),t-^0;and, 

(Vi) / -r„ ~ i^2,.3(/;l/n) , if S. 

Proof: As: 

f(x,y,z)-T„(f;x, y) 

= ^i^^i,^.^„^k(f)^^PWKx + kyy + k^z)}X9i\ k I /n)^'")/(l + 0(|^ I /rtf'") 

+ S;t>„ i^k (/) exp{/(^^x + y + }), 

/-r„(/) < f-g + g-TM + T,{g-f\ 

<2'/-g +max(l,0)n'^'” g 2 ^. 

2 

Taking infimum over g e ^ > we have 

/-r„(/) <Meir2,.3(/;i/n), 

where Mq = max(2,0). As, 0, t — ^ 0, (i) follows. 

If / G L 2 „^'i^ 2 m > previous Lemma, / - T„ | = 0(n n — > oo . 





Conversely, assuming / -r„ = OCn""'") , if / is a fixed natural number, for all n > / 
sufficiently large, there exists a constant M such that 

- Si<;t</ (/)exp{z-(*^x + kyy + k^z)}) 0(1 k 1 1 rtf'" _ I + 0(| A: | ' 

- ^i<k<i(^k(f)^^p{KKx + kyy + k^z))})(d\kf'")l + dik/nf'"]~^ . 

Taking limit as n 0 ° , we have 

i^k (/) exp{i(A:^x + kyy + k^z )) })(0 ) < M . 

Letting / ^ 0 ° , we conclude that / 6 • This completes the proof of (ii). 

For (iii), if / is a constant, / - r„ = 0 = o(n ~^"‘ ) . Conversely, if / -T„ = o(n"^'” ) , 
as in the proof of (ii), we get 

'Lk>M(f)^^P{i(k^x + kyy + k^z)}\dk^"‘) =0. 

Hence, (/) = 0, for all k >1, and / must be a constant. 

To prove (iv)-(vi), if K 2 ^ j(f-,t) = Oi\ir(t^’")),t->0, f -T„^ = 0(i/r(n-^'” )), n 0 = , 
and if isr 2 „, 3 (/;t) = o(v/'(r^'")), t — > 0, '/-r„ | = o(vr(n"^'")), n — > 00 , follows from 
(i). For the rest, assuming f -T^ = 0(iir(n ~^’” )), n — > o® , 

(f - r„ (/);t) + (T, (/ - g);r) + (T„ (g);0 

<0{win^")) + t^"' / -^l+ ^L-- 

Taking infimum over g e L 2 „^ 2 , 2 m^ ’ 

and the results ifollow from Lemma 1.2.2. 

Taking \i/(t) = , for the Bemstein-Lipschitz orders n"“ , we have: 

Corollary 2.4.4. 

Let / £ L 2 „ 3 ^ , and 0<a<2m. Then: 

(i) /-T„ =0(n"“),n-»«>, iff = 




(ii) |/-r. = 0 (n-“),„->=o, iff 

(Hi) f-L -n-^ if K„,,(f-.t)~l‘. 

In the MRI reconstructions using FFT, the relevant choice cx^ = 2~" satisfies 

- (l/2)a„ , so that in the converse parts a given degree of approximation along 

the sequence {2"} itself implies the corresponding smoothness of/ as given in the 
theorem and the corollary. 



2.5. NUMERICAL SIMULATIONS AND DISCUSSION 


In this section, numerical simulation results of approximation by two- 
dimensional Butterworth filter are presented. All the experiments in this section are 
done on test image set II. 

Figures 2. 1-2.3 show segmented images obtained by the application of the 
multi-resolution MLE-based segmented images on PD, T 2 and Ti weighted images 
approximated with two-dimensional Butterworth filter with 0 = 1, 2, 4, 8, 16, 32 for 
the values of m = 1, 2 and 3. The number of iterations taken with 0=1,2, 4, 8, 16, 32 
for m = 1 are 57, 86, 59, 32, 45 and 70 respectively. The number of iterations taken 
with 0=1,2, 4, 8, 16, 32 for m = 2 are 73, 49, 50, 53, 57 and 72 respectively. The 
number of iterations taken with 0=1,2, 4, 8, 16, 32 for m = 3 are 53, 53, 53, 54, 49 
and 50 respectively. 


0 

m=l 

m=2 

m=3 

1 

0.723859 ' 

' 0.091625 

0.015913 

2 

1.314130 

0.177635 

0.031546 

4 

2.249419 

0.335402 

' 0.062014 

8 

3.591323 

0.606762 

0.120002 

16 

5.362906 

1.033300 

0.225807 

32 

7.584896 

1.637128 

0.406243 


Table 2.1: % Relative LI Errors in the Approximation of Original PD weighted Image with 2-D 

Butterworth Filter 


0 

m=l 

m=2 

/n=3 

1 

0.968174 

0.142415 

0.026660 

2 

1.730013 

0.274785 

0.052799 

4 

' 2.890335 

0.514285 

0.103590 

8 

4.462893 

0.916510 

0.199689 

16 

6.385664 

1.525403 

0.373099 

32 

8.597831 

2.343898 

0.662873 


Table 2.2: % Relative LI Errors in the Approximation of Original Ta weighted Image with 2-D 

Butterworth Filter 


e 

s 

II 

m=2 

m=3 

1 

0.719274 

' 0.091599 

0.016963 

2 

1.309699 

0.176995 

0.033582 

4 

2.256373 

0.332299 

^ 0.065835 

8 

3.645170 

0.596168 

" 0.126734 

16 

5.534366 

i.005681 

0.236241 

32 

7.991690 

1.582371 

0.418477 


Table 2.3: % Relative LI Errors in the Approximation of Original Ti weighted Image with 2-D 

Butterworth Filter 


e 

m=l 

m=2 

m=3 

1 

0.841595 

0.089839 

0.014514 

2 

1.558778 

0.175110 

0.028807 

4 

2.749586 

' 0.333892 

0.056763 

8 

4.562928 

0.614222 

0.110335 

16 

7.095091 

1.073159 

0.209350 

32“ 

10370586 

1.759173 

0.382140 


Table 2.4: % Relative L2 Errors in the Approximation of Original PD weighted Image with 2-D 

Butterworth Filter 


9 

m=l 

m=2 

m=3 

1 

0.897747 

0.116682 

0.021314 

2 

1.634145 

0.225710 

0.042227 

4 

2.810198 

0.424509 

0.082908 

8 

4.512003 

0.763304“' 

0.160048 

16 

6.748191 

L289840 

0399839 

32 

9.461606 

2.029814 

0.535309 


Table 2.5: % Relative L2 Errors in the Approximation of Original T 2 weighted Image with 2-D 

Butterworth Filter 


9 

m=l 

II 

s 

m=3 

1 

0.901670 

0.089823 

0.015247 

2 

1.681128 

0.174589 

0.030208 

4 

2.997882 

■ 0.331448 

0.059309 

8 

5.050265 

0.606569 

' 0.114506 

16 

7.985939 

1.057225 

0.214654 

32 

11.858933 

1.744758 

0.384177 


Table 2.6: % Relative L2 Errors in the Approximation of Original Ti weighted Image with 2-D 

Butterworth Filter 



The percentage relative errors in the approximations of the two-dimensional 
Butterworth filter applied to PD, T 2 and Ti weighted images with 6=1, 2, 4, 8, 16, 32 
corresponding to the values of m = 1, 2 and 3 are given in Tables 2.1 - 2.6. 

Graph 2.1 (% Relative LI errors in Approximation with two-dimensional 
Butterworth filter (m=l,2,3)) and Graph 2.2 (% Relative L2 errors in Approximation 
with two-dimensional Butterworth filter (m= 1,2,3)) present the behavior of % relative 
errors with respect to the values of m in each of the weighted images. It is clear from 
the plots that smoothing increases with the increase in the value of 0 and decreases 
with the increase in the value of m as mentioned in section 2.1. 

In all the segmented images, the total number of segments obtained are 
thirteen and visually they are indistinguishable from each other. Out of all these 
segmented images, one would select that particular image which gives the better 
contrast between the segments and the overall sizes of each and every segment 
matches with that represented by given original PD, T 2 and Ti weighted images. 
White matter and gray matter are well separated by a boundary in each of these 
images. 
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Figure 2,2; MLE B 
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CHAPTER III : DISCRETE CONVOLUTIONS FOR 

SEGMENTATION 


This chapter is a study of approximations and the resulting segmentations by 
discrete convolutions with Korovkin, Jackson, generalized Jackson and De La Vallee- 
Poussin filters. Sections 3.1-2 contain the basic results required for the subsequent 
analysis. Section 3.3 contains the direct and inverse results on the approximation pf 
Korovkin, Jackson and generalized Jackson operators. Discrete convolutions with the 
De La Vallee-Poussin kernels have been studied in Section 3.4. 

Numerical results corresponding to Korovkin, Jackson and generalized 
Jackson operators, along with their comparisons, have been considered in Section 3.5. 
Section 3.6 includes numerical simulations with De La Vallee-Poussin kernels. 


3.1. INTRODUCTION 

Convolutions such as with a Gaussian filter have discretization errors in their 
implementations. This necessitates a further error analysis and care during their 
implementations. Also, the Gaussian filters G{X,a)do a spherically symmetric 
isotropic filtering. 

However, the rectangular image size, as well as the nature of noise in the 
phase and the frequency directions in MRI may necessitate applying cartesian 
products L;(L^(L„(/;x);y);z)(3-D), L„(L„(/;x);y) (2-D), and, L„(/;x)(l-D) of 

possibly different filters in the x, y, and z directions. 

The degree of approximation results in such separable cases are essentially 
one dimensional. For example, in 3-D with 

'/ =max{|/(x,y,z)|:-7r <x,y,z<;T}, 


we have: 



xuu 


Lemma 3.1.1. 

If the convolution filters L/ , M„f , NJ -»/, for feC^^, 

m^n (/) ~ f ~ + 0(1/ ^2(^)) + 0(1/ (p^in)), {l,m,n — > i») , 

iff 

~ f\\ = 0(1/ (Pi(l)), (/-»oo), 

=0(l/(P2(m)),(m^oo), 

=0(l/(P3(n)),(«-^oo). 

Proof: Since the approximation property for C 2 „ implies uniform boundedness, let 
the constants L, M, N be such that L, <1, <M , <N . Then, the 1-D 
results imply 

Li (M (N„ (/; x); y); z) - / (x, y, z) = Li (M^ (/V„ y); z) - f(x, y, z) 

< Li(M ^{N^if(^,'n,cy,xy,yy,z)- f(({x,r]Xy,y);z) 

+ Li {M„ if (x,i], C); y); z) - fiix, y, C): z) 

+ IIA (/ (x, y,0\z)-f (x, y, z)|| 

= LM0(1/(P3 (n)) + M 0 (l/(p 2 im)) + 0(1 /(pi (/)) 

= 0(1/ (Pi (/)) + 0(1 /(P2 (m)) + 0(1 /(P3 in )) . 

Conversely, the 1-D results follow from the 3-D relation by letting m, n I, n 
and, /, m — > oo. 

In view of this Lemma, we shall consider a class of discrete 1-D convolutions 
associated with continuous convolutions with trigonometric polynomials that can be 
easily implemented both in the spatial and frequency domains in any of one, two, or 
three dimensions and observe that they possess essentially the same approximation 
and filtering properties as the parent continuous trigonometric versions. Moreover, 
their implementation in the case of two and three variables can be done separably. 

Some continuous trigonometric convolution operators are; 

De La Vallee-Poussin: 

Vnif;x) = [(2n)!!/(2;r(2n-l).^0^J_^_^/(0cos2”^^^^ V 



Korovkin: Fejir-Korovkin: 


, Jit) ^ 

i cos{t-x)-cos{nln) J 


Jackson: 


W2./w(/;^)=j ,/(r)3/(27i>T(2n2+l))j®*‘^K^ dt. 

■[ sin[(r-;c)/2] J 


Generalized Jackson: 




where A„ 


[-7t,7t\-^ I sin[(r-x)/2] J 
[sin(nr / 2) / sin(r / 2)f’’dt. 


''P-P .[-7:.7t] 

For these there hold the asymptotic formulae (see, e.g., P.P. Korovkin [1960], I.P. 
Natanson [1964], R.A. DeVore [1972], R.K.S. Rathore [1974a]): 

V„if;x)-f(x) = f '(x) / n + o(l / n), n -4 oo, (/ (n) > n + 2), 

A -2 iftx)-f{x) = n ^f'(x) /(2n^ ) + o(l/ ), n -> <», (/(«) > n), 

L 2 n -2 (ftx) - fix) = 3f\x)/(2n^) + o(l/«'), n ^ oo, (/(«) > 2n), 

Lnp-p if\x) - fix) = (l - / Po^"'’"^^ )nx) + o(l/n^), n -4 oo, {l{n)>np-p■^r2), 

where [sin(/zt 12)1 sin(t / 2)Y^ = (1 / 2) + Xi<i(:<np-p . 

Another class of useful convolution filters is given by the Poission operators: 

(1-r^) 

^(/;;c) = [l/(2;r)] fit)^- ^ \ \ 

.[-»,«] (l-2rcos(r-,v) + r^)j 

which are positive for 0 < r < 1. 


3.2. DISCRETE CONVOLUTIONS 


For purposes of segmentation it is desirable that the filtering does not lead to 
extraneous, or artifactual segments. This suggests the use of linear positive filters 
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which by definition do not introduce a zero crossing if the function on which they are 
applied are non-negative. For convolutions: 

K (/; Jf) = fit)K^ (t - x)dt, 

with general even trigonometric kernels of a form 

Kn (x) = i ^/2 + cos 0, 

P.P. Korovkin [1960] obtained Pj^ — >1 as the condition for the convergence 
L„f —^fE C 2 J 1 and the following quantitative error estimate: 

Ln(fix)-f(x) < {^ + mn[0.-Pi J/2f^)(o(lJm), 
with m any positive number and (o(S) the modulus of continuity of the function 
fix). 

In our treatment below we restrict to the case of linear positive filters. One of 
the cases, namely that of De La Vallee-Poussin operators that we take up for a 
detailed study actually does much more. I.J. Schoenberg [1959] has shown that these 
operators are variation diminishing, i.e., the number of zero crossings of V„(/) do not 

exceed the number of zero crossings of the function / for any / e C 2 „ ■ The 
convolutions with the Gaussian function G(X,cr) are also known to possess this 
propertry. It is therefore no wonder that these convolutions seem to provide a 
satisfactory usage in image segmentation. 

Let / e C 2 JC , the space of 27r-periodic continuous functions and the ordinary 
modulus of continuity of/ be defined by: 

coif,d) = max |^_^|^5 I fix) - fiy) |, (5 > 0). 

R. Bojanic and O. Shisha [1974] considered the discrete convolutions 
if ’X) [2 /(wi(?z) + y (^A,)7i(n)+2 

associated with the continuous convolutions: K„(/;x) = (I/tt) f{t)^^{x-t)dt, 
with non-negative trigonometric polynomials 



where the nodes of the discrete convolution, in general, are given by 

h,i(n) = -^0 + 27dcfl(n),l <k< l{n). 

R. Bojanic and O. Shisha [1974] considered the case jCo=0, and obtained the 
quantitative error estimate 

^n.m(«)+2 (/; ^) - fix) = (1 + 7r)Co{f-,[il - Pj „ )/2f^). 

In the following we shall study the more general discretizations 

^n,l(.n) ifyX) ~ -j:), 

with ;co as any fixed point, and l(n)>m(n) + 2 a positive integer. The basic idea 

behind these discretizations is the trigonometric polynomial precision ‘n’ of a 
Gaussian quadrature based on ‘n+l’ equidistant nodes: 

.l-n,n] ^ )> i^k =XQ+2nk /{n + l),0<k< n), 

where Xg is an arbitrary point and (x) is any trigonometric polynomial of order n. 
Accordingly, if /(x) is a trigonometric polynomial of order p, 

Kn.nn)ifyX) = ^nifyX), 

provided /(n) > m(n) + p. Using it, a generalization of the R. Bojanic and O. Shisha 
[1974] tu(/;<5) error estimate, and also an coif '-,6) error estimate is: 

Lemma 3,2.1. 

Let /(n) > min) + 2 . Then: 

(i) if-,x)- fix) < (1 + 7t)(0if-,[il - Pi,„ ) / 2f^ ), / e ; and 

(ii) K„^n,^if-,x)-fix) <7t^ /'(xX(l-Pi,„)/2 + 2;r2[(l-Pi.J/2]''2 

x£u(r;[(l-p,,J/2f2),/'eC2^. 

Proof: As the operators are linear positive (R.A. DeVore [1972 Th. 2.4]), for 

/ ^ ^271 1 

X tu(f; (sin^ [(r - x) / 2]; x)), 

and,for /6 Cj;,, 
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Knm^f\x)-f{x) < fix) l-K,j^„^il-,x) + fXx) (K,j^^^i&mit-xy,x) 

+ (sin^[(r - ^)/2];;c))+;r2 (l + ( 1 ;^ 

X fc./(n) (sin^ -x)/2];x)y\Q)it (sin^ [(r -x)/2]-, x)). 

Since (? — x), sint(j)^(t — x), cost (j)^it — x) are trigonometric polynomials of 
orders min), m(n) + l, and, m(n) + l, respectively, for /(n)>w(n) + 2, we have: 

~ (^tn t, jc) — sinj:, K^ i^^-^icost',x) = Pi^n cosjj. 

Hence, 

K , x) COS X isr^j^^^,j^(cosrix)sin jc = 0, and 

^n,i{n) (sin ^ [(t - Jc) / 2] ; x) = (1 / 2)[1 - (cos t; x) cos x - (sin f, x) sin x] 

= (l-Pi.J/2, 

from which the Lemma follows. 

Also, one has the following asymptotic error formula in the approximation of 
fix) by the discretizations K„n„^if;x): 

Theorem 3.2.2. 

If lin) > min) + 2, /is a bounded 2;r-periodic function, and f" exists at a point x, 
Knm U\X)- fix) = (1 - Pi,„ )/'(x) + 0(1 - p,,„ ), 
iff for some k>2, lim„^„ (1 - p^t, J/(l - Pi,„ ) = 

Proof: Both sin 2t <j)„ it - x) and cos2r^„(r-x) are trigonometric polynomials of 
order min) + 2. Hence, if /(«) > min) + 2, also there hold: 

(sin 2r; x) = P 2 .„ sin 2x, (cos 2t\ x) = P 2 .„ cos 2x. 

From this the result follows since for a sequence L„ of linear positive operators there 

holds (R.K.S. Rathore [1974a, Th. 1.2. 1.5]) the asymptotic relation 

L„ (fix) - fix) = fix){L, (l;x) -1} + f'ix)L„ (sin(r - x);x) 

+ [2/'(j:) + o(l)]L„ (sin ^ (r - x) / 2; x), 

as n — > oo, iff for some k = 2,3, 

lim„_^„ (L„ (sin ^ kit- x) 1 2; x))/ iL„ (sin ^it-x)/2;x)) = k^. 
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It may also be noted that the asymptotic formula holds uniformly in x on each 
closed sub-interval of an open interval on which /'exists and is continuous. In 
particular, for twice continuously differentiable functions in € 2 ^ the asymptotic 
formula holds uniformly in x on the real line (R.K.S. Rathore [1974a, Th. L2.1.5]). 


3.3. DEGREE OF APPROXIMATION BY DISCRETE 
CONVOLUTIONS 


In particular, for the De La Vallee-Poussin, Korovkin, Jackson, and the 
generalized Jackson continuous convolutions given above, their discretizations are 
given as follows: 


‘ 2 


'-‘*,/(n))COS (l 


kyl{n) 




( sin(;r/n)cosn [(?4 ;(„) -x)l2] 
\<.k^l{n) k,Kn) ^ COSff^^ - x) - COSf^T / n) 


x)/2,, 
2 ~ 


^2n~2,Kn) (/ 5 X) = [2//(n)]Ei^,^;(„)/(?*,/(n))[3/(2n(2n2 +1))] 


sin[n(t^_y(„) -x)/2]. 
\ sin[(ri_,(„) -x)/2] ^ 




;sinW<.,„.,-j:)/2];/ 
sin[a,,,(„) -x)/2] ^ 


For these discrete filters, there hold the asymptotic formulae: 

(f‘,x) -f(x) = f\x)ln + 0(1/ n), n (/(n) > « + 1) , 

K-2,i{n) fix) = ITT V'W/(2n^) + 0(1 /n^), n (/(n) > n ) , 

^n- 2 ,iM ~ '^f"(x:)K‘2n^) + o(l/ n^), n — > Qin) > 2n ) , 

f Q (.np-p) > 

Vp,/(n(/;^)-/W= \np-p) f(x) + 0 (l/n\n^oo,(l(n)>np-p + 2).- 

I Po J 

These asymptotic formulae are the same as those for the original operators 
(see e.g., P.P. Korovkin [1960], I.P. Natanson [1964], R.K.S. Rathore [1974a]). To 



study the direct and inverse results on their degree of approximation, some definitions 
needed are as follows: Let as usual denote the space of 2;?r-periodic continuous 

functions on R, and C 2 ^ the sub-space of functions having continuous second order 
derivative. Let / = max/(x) , /£ C 2 „, and let |4 = g| + V|, geC^^^Lttthc 
Peetre’s K-functional for / £ be defined as: 

/:(/;?) = inf{|/-g +t^ g'\\: 8 GC 2 „^},i> 0 . 


Lemma 3.3.1. 

If l(n) > m(n) + 2, for some constant M, 

S ^n,l{n)S — (1 ~ Pl,/i ) 2> 8 ^ ^2jt • 

Proof: By a)(/'; 5) error estimate, 

Knjin)(8>x)~g(x) <7r^g'(x)(l-p,„)/2 + 2;r^[(l-Pi„)/2]‘^2 

X£0(g';[(l-pi„)/2]''2)<((l-Pi,J/2)(;r2+2;r2)lg'||, 

from which the result follows. 

Lemma 3.3.2. 

If l(n) > m(n) + 2, 

(8>^) < (7an(n)/2f(l- pi J g"' 2 , g e . 

Proof: Choose a such that g'(o:) = 0. Then, 

■ (s;:t) = «(a) + (« - x)^ gW'2;x), 

SO that by the Bernstein’s inequality, 

<[m(n)f ii:„,,(„)((t-x)2g'(|)/2;x) 

< [mn(n)/2]^ (2sin^ (r - x) / 2; x) \\g\ 

<[mn(n)/2f(l-Pi, Jig'll- 



Theorem 3.3-3. 

Let f{x)eC 2 „, and xifeS. If =0([m(n)r2) ~ l/„2 l{n)>m{n) + l, 

then: 

0) ^n,i{n)f~f - — > 0, n — > oo, where A is a constant; 

(ii) ^«,/(n)/ - / = Oiy/in'^)), n-^oo, iff K(f-,t) = 0 ; 

(iii) ^n,i(.n)f~f =o(y/(n~^)),n-^oo, iff K(f\t) = 0{y/(t^))j -^0-, 

(iv) ~ K(f-,n~^), if K{f-/’^)e S. 

Proof: For (i) we have 

^n,l(.n)f~f - 5“/ + ^n,l(.n)S~8 + ^n,H,n)if ~ S) 

<2f-g +M(l-p,,„)g.2. 

Taking infimum over g e € 2 ^'^, 

where A = max{2,Af }. 

To prove (ii)-(iv), using (i), if ^ 0, 

Kn.l(n)f-f =0(l/r(l-Pi_J),n^oo, 

and, if, K(f-,t) = o(y/it^)),t-^0, 

Knmf-f =oW-pi.J),n-^~. 

For the inverse parts of (ii)- (iv), using the previous two Lemmas, 

K{f-,t) < K(f - A:„./(„)(/:0) + K(Knjin)(f - + 

— f ~ ^n,l(n)f I •^n,/(/i) (f ~ 8) ^nJM 8 | 

< +t^i7an(n)/2)^{ f - g +[(1-Pi,„)]^ 2 }’ 

2 

Taking infimum over g e C 2 „ , 

and the results follow from Lemma 1.2.2. 
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For the familiar Bernsteinian orders n “ , 0 < a < 2, we have the immediate: 

Corollary 3.3.4. 

Let f(x)e 0<a<2, 1-Pi„ =0([m(n)]~^) ~lJn^ and l(n)>m(n) + 2. Then: 

(i) , ^„./(„)/ - / = O0r“), n -» oo, iff K(f-t) = 0(r“),r 0 ; 

(ii) =o(n““),n-»oo, iff /,:(/;;) = o(r“Xr-> 0 ; 

(iii) if ~ n'“. 

For the De La Vallde-Poussin kernel there holds: 1- pj^ = 0(l/n). Hence the 
above results are not applicable to the discrete filters which we shall 

take up separately in the next section for a much more detailed study. 

However, since the relation 1-Pi,„ = 0([m(n)]~^) ~ 1/n^ is valid for the 

earlier mentioned Korovkin, the generalized Jackson, and Jackson kernels (the case p 
= 2, in the following), for the discrete convolutions L„p_pi^^j(f;x) 

and L 2 „_ 2 j(n) (/;^)> we obtain the following: 

Corollary 3.3.5. 

Let f(x) E C 2 „ , and y/e S. Then, if l(n) > n : 

(i) \- 2 ,i(,n)f ~ f ^ A(O2if',l/n)-^0,n—^°°, for a constant A; 

(ii) =0()/r(n-2)),„->oo, iffru2(/;t) = O(vr(r2)),?_>0; 

(iii) A„_2,;(„)/-/ =o(V^(n“^)),n-^oo, iff (y2(/;r) = o()/r(r^)),r->0; 

(iv) ^n-2,l{n)f ~ f ~'^2^f'’^ )> if ^2(/>^ ) S 5. 

Corollary 3.3.6, 

Let fix)E C 2 „, p>2, and y/e S. Then, if lin)>np- p + 2: 

(i) Vp,/(„)/-/ ^^ 0 ) 2 (/;!/«)-> 0,n->oo, for aconstantA; 

(ii) VpA«)/-/ =0(V^(n'" )),«->“, iff co2(f-,i) = O(y/{t^)),t^0-, 

fe '.i 



(iii) ^np~p,i(n)f f — o(i//'(n = 

(iv) Lnp-pmf-f ~ if (0^{f-t^'^)ES. 

In particular for the Bernsteinian orders n"“ , 0 < a < 2, there hold: 

Corollary 3.3.7. 

Let fix) € , 0 < a < 2, and l(n) > n. Then: 

(i) =0(n““),n-^oo, iff (02(/;t) = O(r“),r->0; 

(ii) \-2,i(n)f~f =o(n “),n->oo, iff (D2(/;r) = o(t“),r— >0; 

(iii) V 2 , /(,!)/-/ ~ n"“,if ft) 2 (/; 0 ''t“. 

Corollary 3.3.8. 

Let / (x) e C 2 fc , 0 < a < 2, p > 2, and l(n) >np-p + 2. Then: 

(i) ^np-p,i(.n)f ~ f = 0(^ ^ iff £02(/;0 = 0(r'^), r — > 0 ; 

(ii) ^np-p,i(n)f - f =o(n‘“),n-^oo, iff co 2 (/;t) = o(r“),t-»0; 

(iii) ^np-p,l(n)f ~ f ~ ^ if <^ 2(/>0 ~ 

In these corollaries, we have used the well-known equivalence 

AKif-,t) < (02if;t) < BKif-t), it > 0) 

of Kif\t) and 0)2 if',t ) , where A and B are positive constants, and indeed the second 
order modulus of smoothness is defined as 

(O 2 if\t)= supi^i^, (sup^[_^_^j I /(x + /i) - 2/(x) -f- /(x - /i) 1 ), / e € 2 ^ . 

3.4. DISCRETE CONVOLUTIONS WITH DE LA VALLEE- 
POUSSIN KERNELS 

In this section we study discrete convolutions with the De La Vallee-Poussin 
kernels, whose continuous version is given by: 



V, (X) = V, (/; 4 = ‘/(Ocos^”'-' 

27i(2n-l)\\:^ ^ 2 j 

Some well known results (I.P. Natanson [1964]) about the continuous convolutions 
V„(f-,x) are: 


(i) If a function / (x) e has modulus of continuity Q){6) , then the inequality 

K, W - fix) <3(0 

^■slnj 

holds for all values of x. 

(ii) If = sup{max'l/,(x) — /(x)|}, wherein /(x) runs through all functions 

with period In of class Lipitt , then 


where lim p„ = 0 . 

n-^oo 


UniCC) 


2“ p„ 

^Jnn°‘ K ^ 


(iii) If the function / (x) 6 € 2 ^ possesses for a value of x a finite second derivative 
fix ) , then the formula 

n n 

with lim p„ = 0 holds true for this x. 


(iv) If at a point x , fix) 6 ^2n has a finite derivative f'{x) , then the relation 

lim !/„ (x) = f{x) 

AZ— >00 

holds true for this value of x. 

Generalizations of some of the above results by P.C. Sikkema and R.K.S. Rathore 
[1976] are: 

(v) Let 0 < a < 1 , and p be a positive integer. If 

f~ p '' 

E((x,p-,V„,x) = sup< V,„ fit)-X J i^-^y >, xeR, 

f . v=o V. _ ^ ^ 



where / runs through all functions with 6 Lip, a, then with 
{a)p =nf=i(a + /-l), 


2P+2a-l jp^ p+a+l ^ 

"Tl + a), 


o ) £1“ 9P+ap/p+a+K 

/i- < lim« ^ E{a,p-,V„,x) < ^ . 


(vi) If p is a positive integer, then 


lim n ^ £ (1, p ;y 

n oo ^ n 


x) 


2 + 1) 
( p + 1)! -JiT 


(vii) Let m be a positive integer. If /e Q^and at a point ;c6R, f'^\x) exists, 
then 


(viii) y„(sin^*r/2;0) 
3.2.1.]) 


n-^cix 


= 7r n *,/:> 1. (V.S.N. Kaliprasad [2000 Lemma 

V 2 j 


Lemma 3.4.1. 

Let m > 0, n > 1 be positive integers, and D = d/dx. Then 

where [m/2] is the integral part of mil, and P^jin ) , (0 <y < [m/2]), is an algebraic 
polynomial in n of degree < j. 

Proof: For m = 0, the result is trivial, with q (n) = 1 • Since, 

d(cos ^‘{x/2))= (cos^" (xl 2)\-n) tan(x/ 2), 
and 

L»2[cos2"(x/2)] = [cos^”(x/2)][n(n-l/2)tan2(x/2)-2n], 

the result is true for m = 1,2. Hence assuming it for m, 

D'""* (cos^'’ (X/ 2)) = feo, (n)D^os^" (xl2) tan'”'^^’ (x/ 2)}) 

= cos"” (x/2)] 2o,;<(m/2] («) /2-j)-n) tan'”"'-"^' (x/ 2) 

+ (m/2-7)tan"'"‘-"^^'"'^(x/2). 



= (cos 

+ (m/2 - ;■ + .,1 (n)]tan''‘"'-2^(;c/2)} , 

and Pm+\^m-¥\-j (”) = [Pm,m-j (n){{ml2 - j) - n) + (ml 2- j + 1)P„ ) is a polynomial 

of degree < m+l-j. Hence the result holds for m+1, completing the proof. 


Lemma 3.4.2. 

There exist trigonometric polynomials (x) , 0 < yt < 2p, of order p such that 

(■^) = ) ,0 <k< 2p, p > 1 . 

Proof: Let be the functions defined by (R.K.S. Rathore [1974b]) 

where \t,x), k = 0, 1, ... , 2p, are the determinants of order (2p+l) given by 
(t,x) = (a,y , where with [//2] denoting the largest integer not greater than j/2, 


and 


a[f = ^ sin[;72]L 
cos[;72]r, 


7=1 

j is even , 
j is odd 


■*•5 

= < sin[;72]x, 
cos[772]x, 


7 = 1 

j is even, 
j is odd 


k^O , 


d 


/-I 


- ^i-x ^ifc+l y , 1 ^ i ^ + 1 ■ 


Since the functions Q!'’^\t,x) satisfy 


(t, x) = (t-xf +0(\t-x ) , 


taking 7]^ (x) = (^) = (^>0) > ibe result follows. 


In conjunction with Taylor’s expansion, as an immediate consequence of the 
Lemma, we have the following; 
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Corollary 3A.3. 

If at a point ‘x’, (x) exists where /(t) is a bounded 2;r-periodic function, then: 

fit) = / (x) + (/ ^ ix)/k\)T, it-x) + h,it)it-xy\ l<m<2p, 

where it) , with ix) 0 , is bounded and is continuous at r = x. Moreover, if 
f^ ^ ix) exists at all x and belongs to Cj® for an any e > 0, there exists a 5 > 0, 
independent of x, such that h^it) <e ,0< |r-x| < <5. 

The next result establishes the basic approximation of pointwise derivatives by 
the derivatives of the discrete De La Vallee-Poussin operators V„(/;x) : 

Theorem 3-4.5. (Approximation of derivatives). 

If / is a bounded 27r-periodic function such that for some natural number m, and at 
some point jc, f^"'\x) exists, then with l(n) > n+[ml2]+\: 

Moreover, if (x) exists for all x and belongs to the convergence is uniform 
in X. 

jProof: Using the Corollary above, where we choose 2p = m, if m is even, and 2p = 
m+1, if m is odd, we have: 

Since T/s are fixed trigonometric polynomials satisfying 7)'" (0) = m\5 j^,0<j<m< 
2p, where Sj,^ is a Kronecker delta, 


Hence, it remains to show that the contribution of the o-term tends to zero. For this, in 
view of uhe Corollary above, given an arbitrary e > 0, we can find a <5 > 0, such that, 

for some constant A, and for all . Hence, with = m, if m is even, and fi = m+1, 

if m is odd, in view of the 2;7:-periodicity of f (/) using which we could assume that 
(x - ^ we have, for /(n) > n-pi/l+m+l = n+[m/2]+l, 

[(2n)!!/(/(n)(2n-l)!!)E.,,,,,,o(x-t^ 

^ Sos2 j<.{ mi2\ 00 tan ((.V - f f. ) / 2) 

“ ^^022 ji[m/2] (”) ^n-ii/2+j,l(n) Sin ((.JC - ) / 2) 

+ (A/5^) ((x - )/2y,x)) 

by (\'iii) above, from which, due to the arbitrariness of £ > 0, it follows that the 
contribution of the o-term tends to zero as « — > <= 0 . 

The uniformity of the convergence: 

in the case when f^'"^ (x) exists for all x and belongs to Cm, is clear since then A and 
the S in the above could be chosen to be independent of x. 

The next result provides an asymptotic formula of Voronovskaya type in the 
approximation of a derivative at a point: 

Theorem 3.4.6. (Voronovskaya Type Asymptotic Formula). 

If/ is a bounded 2;r-periodic function such that for a natural number m, at some point 
X, (x) exists, then with l(n) > n+[m/2]+2: 

lim„^ - /‘•>(2:))=/‘"+«W. 

Moreover, if (x) exists for all x and belongs to € 2 , 1 , the convergence is uniform 


m X 



Proof: Using the Corollary above, where we choose p = [(m+3)/2], i.e., Ip = m+1, if 
m is even, and 2p = m+3, if m is odd, we have: 

+ 0(X - ) X COS^" (X - ) / 2 

^ ^0<2 y<[m/2] (”) ■' ((x - 4,/(n) ) / 2)_ 

= Ei<;^m +2 0’!)~7^'^ (^y (^ -^);^) +[(2n)!!/(/(n)(2n - 1)!!)] 

^ ^l<kil(n) pOS " ) /2_ 

^ ^ 0 S 2 y<[m/ 2 ] ^m,m-j M tan ((x - ) / 2)_. 

First of all we consider the contribution of the o-term, for which given an arbitrary e> 
0, we can find a <5 > 0, such that. 


for some constant B, and for all . Hence, with p = 2[(m+l)/2], i.e., ju = m, if m is 
even, and p = m+1, if m is odd , in view of the 2:?r-periodicity of f(t) we could assume 
that (x-r^_;(„)) <7r, so that x-t^n^^) <7rsin((x-r^j(„))/2) . Hence, for l{n) > n- 


\x/2+m+2 = n+[m/2]+2: 


[(2nV.\mn)(.2n - 1) oU - (, cos^” 


X 


=o(S, 


^0<2y<[m/2] 




P . (yj) 

0<2 j^[m/ 2] ^rn,m-‘jV‘'J 


X K-n / 2+y./(n) (« sin ((x- ) / 2) + (5 / 5 ^ ) sin ((x - ) / 2); x)) 

by (viii) above, from which, due to the arbitrariness of £ > 0, it follows that the 
contribution of the o-term is o(n^), as n -+ oo. 

Next, using [.] for the integral part, since T/s are fixed polynomials of order p = 
[(m+3)/2] and p+n = [(m+3)/2]+n < «+[m/2]+2 < l(n), 

VnAn) (Tj (t - xy,x) = y„ (Tj it - x);x) = U„ (Tj ;0) . 



Let Tj{x) + cosfac + Zj^sinb:}. In view of the Voronovskaya 

asymptotic formula for the operators V„, 

(Tj ; X) = ao / 2 + Pk,n W cos kx + b;^ sin fcc] 

= ^o/2 + 2^j^^^^(l-ii:^/n + o(l/n))[a^ coskx + b^ sinfcc], 

so that 




("») , 


nJW (Ty.x) = vJ"\Tj(l-xy.x) 

~ Po ^skz + b^ sinfejj 

= r/”> (0) - ln)W cosfe + b, sin fa]}'""’ 

= (0) + (1 / (0) + 0(1/ n) . 


z=0 


m) 

z=0 
+ 0 ( 1 //;) 


Hence, since Tj' (0) -m\Sj,^,0<j<m< 2p, where 6 is a Kronecker delta, we 
have 

= S.sism,2(/-)‘‘/‘"'’W(r/'”’(0)+a/n)r/”*"’(0))+oa/n) 

= / (:c) + (1 / «)/ ("’■'2) ) + 0(1 / n), 

from which the required asymptotic relation 

follows. Noting that if f^"''*'^\x)E C 2 ® the 5 in the above could be chosen 
independently of x and, furthermore, that the magnitudes of all the derivatives 
f^^\x), 1 <y < m+2, possess a common upper bound, the uniformity part follows, 
completing the proof of the Theorem. 


Let C 2 /’, (m > 0), denote the space of 2;r-periodic functions having continuous 
m-th order derivative on the real line. Let ||/|1 = max ^{x}\,f £ Cin- Let a Peetre’s K- 
functional, related with the pair {€ 2 ,!^ ,€ 2 ^*^), for/ £ C 2 /’ , be defined as: 

/C„«(/;') = inf I /<■">-«<"> +<" 


In order to handle the approximation of functions as well as their derivatives 
simultaneously, we would naturally identify /<«> (;,) with /(x) . and, Cr.» with Cr» 


Lemma 3.4.7. 

If l(ji) > n+[w/2]+2, 


where M is a constant. 

Proof: By mean value theorem, for some ^ lying between t and x, 

g{t) = 8 {x) + Sisy<,„+i fe (^) / ;! ){t-x) j +(g ("’+2) ^ 2) !)(^ _ xy^+2 ^ I ^ ^ 


(hi) , 


In the rest of the proof, using 2;r-periodicity of g and we restrict ;c to belong to 

[-n, 7tl and shall regard the functions (f-;r)2 defined by the corresponding value 
(r - x) for I r - X |< ;r . Using p = [(m+3)/2], since 

Tj{x) = xUO{x^P^^), 0<j<2p, p>\, 

we have 

{t - x) 2 =T j(t -x) + 0((t- x) ) 

=T j(t - x) + o(sin^^ ((r - x) /2)) 

=T j{t - x) + o(sin'"'^^ ((r - x)/ 2)), 

where the 0-terms holds uniformly in |r-x| < n, so that in view of \\h'\\ < 7%'% h s 
C2;:^ 

S«) = iW + y;^.^_^,,^“W/;!.rj((-x)+4<”«>o(sin"“((t-x)/2)),|t-x|<s. 

Using the estimates: V'.,,(„,'“'(T,(t-x);x) = ;!5^.„ +o(l/n), 0< ;<m + l, 

in which the o(l/n) term holds uniformly in x, of the previous Theorem, with p = 

2[(m-H)/2], 

= + j<”*» o(sm“«((t-x)/2));x) 



J. XU 


= o(n-') + V„,„„/”'(o(sin“^((,-;t)/2)u)) 

X (cos^" U („) tan (x - tt,„„ ) /2 ' 

<n^M , 

for some constant M, from which the result follows after taking supremum overx 

Lemma 3.4.8. 

If /(n) > n+[<gr/2]+l, g <K^ g^^^ - 5 ^ ^ 2 ^^, (i? ^ 0), being a constant. 

Proof. Let /i = 2[((y+l)/2], and l(n) > n+[^/2]+l. For Tfs with p = /V2, for q>0, using 
Vnj(n)^^^ (1; jj) = 1, and V (Tj (t-xy,x) = 0(1), uniformly in x, we have 

= ^/i,/(«)^^^(AU) + 2i<y2^_j(j!)'^A^''^UFy(t-x)+ g^^^ o(sin^((r-jr)/ 2 ));;c) 

= g^^^ 0(l) + F„,,(„/^^(o(sin?((r-;c)/2)):;c) 

= 1^0(1)+ ((2n)!!/(/(n)(2n-l)!!))5;,..^„,0(sin«((;c-(.,,(.,)/2)) 

><(“s'"(->:-'*2w)/2)Sc«2M,,2|V2W“’'''(*-'«mV2 ^ 
= g <’> (Od) + (( 2 : )/ 2));2r)) 

sKgg^‘>\, 

for some constant Kq. Hence the result, it being trivial for ^ = 0. 

Lemma 3.4.9. 

If l{n) > n+[m/2]+2, < nL , h e C 2 n”'\ m > 0, L being a constant. 

Proof: For l{n) > n+[m/2]+2, as in the proof of the previous Lemma, 


JL X Ji' 


= + />« 0(sin"((/-;t)/2));x) 

= *<”> o(n-') + V'„,„„<"*«(o(siii"(((-4/2))i) 

= K2n)!l/(;(„)(2„-l)!!)]5;.^^^__^o(sta”to-<,,„0)/2)) 

x(cOS (j; “^2 J ((0 )^2 

^Ki~p/2-i+;,/(n) (sin 

^riLh^’"^ , 

where L is a constant. The result follows. 

The following is the main direct and inverse result about approximation of/by 
V (case m = 0) and its continuous derivatives by the m-th derivative 

Theorem 3.4 J.O. 

Let f(x) 6 C 2 /', m>0, l(n) > n+{ml2\+2, and y/e S. Then: 

(i) ^ “> 0. n ^ o=, A being a constant; 

(i>) =0(r(n-‘)),/i^~, iff i:„,„«(/;l) = O«t")),l^0; 

(iii) =o(*2(n-‘)),,.-^~, iff X„.2(/;t) = o(*2(,^)),(-^0; 

0'’) n., 

Proof: (i) 

l'.,, «"/-/'"> S V,, „«(/-«) + V„2„o‘"V-s“ + 

< (1 + - g + n"‘M g . 

Taking infimum over g 6 Czn^'^, 


where A = max { l+K, M). To prove (ii)-(iv), using (i), if 

= =0(rtn-‘)).„-^», 

and, if 

^m,m+2 (/’O — 0(V^(? ^ 0 1 i^n,/(n)^ = 0(l//'(n~^ )), n — > oo. 

For the inverse parts of (ii)- (iv), using the previous two Lemmas, 

s /'"> -v.j,.rf +'H } ' 

Taking infimum over g e for a constant M, 

s m{ /<"> -V., ,„/”>/ + 
and the results follow from Lemma 1.2.2. 

The following consequence of a general interest is immediate: 

Corollary 3.4.11. 

Let f(x) E C2n\ m > 0, 0 < a< 2, and l(n) > n+[m/2]+2. Then: 

(i) =0(n-"'^),n-»oo,iff4:„„,^2(/;,) = 0«“),(-»0; 

(>i) =o(n-“'^),n^oo,iff = 

(iii) V,, - /<"> ~ n-“'^ if ~ 

Finally we relate the AT-functional second order modulus 

of continuity tU 2 (/^'"^;< 5 ) of the derivative by showing that they are 
equivalent. (The case m = 0 of the result is, of course, well known. However our proof 
is valid for all m > 0). 

Theorem 3.4.12. 

For m > 0, and fe C 2 /, 

<Uj (/<">;() £ (/<">;«, (>0. 


Proof; Let g e Then, 

• 0 •'0 

Hence, 


Taking infimum over g € we get 

Conversely, with 5h denoting the central difference operator with step h, and choosing 

g(x) = (2h)~^ f{x + r + s)drds, 

_ // 

W = (2/i)-2 ‘\f>’'\x^r + s)drds 

_ / 

= (2h)~^ d2hf^"‘Hx + s)ds 

• —/l 

It follows that 

and taking h = r/2, we have 

^(,+2) 

Also, we have 

/ (x) - g (x) = / (x) - (2/i)-^ f^"'\x + r + s)drds 

• “h * —h 


= / (x) - (2/z) (x + r + 5) + / (x - r + J) 

^ f(.'n)^x + r-s) + {x-r- s)]drds 

from which it follows that 

4 < (W2(/^'”^;2/z)+<W2(/^'"^;/i).<5t02(/^'”^/i). 


Hence, 


4a:„,„«(/;0S4{ /<”>-«« 


+»" j<9W2(/‘”;*), 


completing the proof. 



As an immediate consequence, the earlier direct and inverse results for the 
discrete De La Vallee-Poussin filters could be alternately re-stated as: 

Corollary 3.4.13. 

Let/fxj e Cin”, m>0,l(n)> n+[m/2]+2, y/ e S, mdO < a<2. Then: 

0) Kpi,/(h) ^ - ■^‘^2 (/>^ *^^) ~> 0, n oo, being a constant; 

(ii) ^n,l{,n) ^ f - =0(i//'(n”‘)), n-^oo, iff ; 

(i'i) =o(i(^(n-')),n^<», iff tti2(/;() = o(rtl')),l^0; 

(iv) V,, = 0 (n-«'^), iffffl,(/;() = 0 ((“); 

(V) =o(„-“'^), iff 

(Vi) ~"-“'^if<B2(/lI)-(». 



3.5. NUMERICAL SIMULATIONS AND DISCUSSION FOR 
KOROVKIN, JACKSON AND GENERALIZED JACKSON 
OPERATORS 


Numerical simulation results for Korovkin, Jackson and Generalized Jackson 
operators with n = 128, 64, 32, 16 and 8 are presented in this section. The experiments 
were done on test image set I. 

In Figure 3.1 (Original T 2 weighted and corresponding images approximated 
with Korovkin operator), the first image is original T 2 weighted image and other 
images correspond to the Korovkin operator with n = 128, 64, 32, 16 and 8 
respectively. Figure 3.2 (MLE Segmented images of original and corresponding 
images approximated with Korovkin operator) contains MLE-based segmented 
images associated with the images in Figure 3.1. The number of iterations taken for 
original image is 59400 and for n = 128, 64, 32, 16, 8 are 10250, 36150, 2950, 2750, 
3000 respectively. 

In Figure 3.3 (Original T 2 weighted and corresponding images approximated 
with Jackson operator), the first image is original T 2 weighted image and other images 
correspond to the Jackson operator with n = 128, 64, 32, 16 and 8 respectively. Figure 
3.4. (MLE Segmented images of original and corresponding images approximated 
with Jackson operator) contains MLE-based segmented images associated with the 
images in Figure 3.3. The number of iterations taken for original image is 59400 and 
for n = 128, 64,' 32, 16, 8 are 14400, 10200, 25900, 9000, 2550 respectively. 

In Figure 3.5 (Original T 2 weighted and corresponding images approximated 
with Generalized Jackson operator), the first image is original T 2 weighted image and 
other images correspond to the Generalized Jackson operator with n = 128, 64, 32, 16 
and 8 respectively. Figure 3.6. (MLE Segmented images of original and 
corresponding images approximated with Generalized Jackson operator) contains 
MLE-based segmented images associated with the images in Figure 3.5. The number 
of iterations taken for original image is 59400 and for n = 128, 64, 32, 16 and 8 are 
8450, 16150, 22150, 9200 and 2900 respectively. 


of iterations taken for original image is 59400 and for n = 128, 64, 32, 16, 8 are 8450, 
16150, 22150, 9200, 2900 respectively. 

The percentage relative errors in the approximation with Korovkin, Jackson 
and Generalized Jackson operators are given in Tables 3.1 and 3.2. 


n 

Korovkin 

Jackson 

Generalized Jackson 

8 

18.385885 

18.120655 

17.940020 

16 

18.066454 

17.019840 

16.342531 

' 32 

15.425109 

13.358624 

11.905014 

64 

11.264680 

8.087324 

6.817662 

128 

4.510457 

3.271766 

1.597230 

Table 3.1: % Relative LI Errors in the Approximation of the Original Cross-Section with Korovkin, 

Jackson and Generalized Jackson Operators 

n 

Korovkin 

Jackson 

Generalized Jackson 

8 

20.188421 

19.894699 

19.693750 

16 

19.834602 

" 18.656040 

17.878233 

32 

16.788383 

14.335914 

12.556752 

64 

11.735953 

7.875902 

6.418863 

128 

4802420 

2.915334 

1.400528 

Table 3.2: % Relative L2 Errors in the Approximation of the Original Cross-Section with Korovkin, 

Jackson and Generalized Jackson Operators 

The case n = 

128 gives an image visually similar to the original and the case n 

= 64 gives necessary smoothing without distorting major regions with all operator 
mentioned above. Blurring in the images increases as n decreases from 128 to 8 and it 
is high for n = 32 downwards. The blurring with korovkin operator is more followed 


by Jackson operator and then Generalized Jackson operator. This is clear from Tables 
3.1 and 3.2 and is also observed from the images shown in Figures 3.1, 3.3 and 3.5. 

The segmented image for n = 128 is similar to that obtained from original 
cross-section, result for n = 64 seems to be satisfactory segmentation of major tissues 



along with the pathology. For n = 32, 16 and 8, the boundary between major segments 
like cranial part of the brain, white matter and edema, thickens and smaller regions 
merge with neighboring larger segments. This is observed that the segmented image 
obtained for n = 16 (Figure 3.2) with korovkin operator matches with that obtained for 
n = 8 with Jackson (Figure 3.4) and Generalized Jackson (Figure 3.6) operators. 

With korovkin operator, segmented image obtained for n = 64, a thin boundary 
encloses the region comprising of edema completely which could is useful for 
medical purposes. Also, for n = 64, 32 and 16, the segmented images show three 
regions corresponding to edema which are represented as different growth stages of 
an edema within the brain parenchyma. Similar kind of splitting of edema has also 
been observed in the cases of Jackson and Generalized Jackson operators for n = 32 
16 and 8. In cases of these two operators, boundary enclosing the edema is observed 
for n = 32. Plots of percentage relative errors in the approximation with respect to the 
values of n are given in Graph 3.1. 




Figure 3.1: Original Tj weighted and Conresponding Images Approximated with Korovkin Operator 










Figure 3.2: MLE Segmented Images of Original and Corresponding Images in Figure 3.1 










Figure 3.3: Original Tj weighted and Corresponding Images Approximated with Jackson Operator 










Figure 3.4: MLE Segmented Images of Original and Corresponding Images in Figure 3.3 









Figure 3.5: Original T 2 weighted and Corresponding Images Approximated with Generalized Jackson Operator 









n=16 


n=8 


Figure 3.6: MLE Segmented Images of Original and Corresponding Images in Figure 3.5 
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Graph 3*1; Plots of % Relative LI and L2 Errors with respect to the Values of n 



3.6. NUMERICAL SIMULATIONS AND DISCUSSION 
CORRESPONDING TO VALLEE-POUSSIN KERNELS 


Numerical simulation results for De La Vallee-Poussin operator with n = 
4096, 2048, 1024, 512 and 256 are presented in this section. The experiments were 
done on test image set I. 

In Figure 3.7 (Original T 2 weighted and corresponding images approximated 
with De La Vallde-Poussin operator), the first image is original T 2 weighted image 
and other images correspond to the De La Vallee-Poussin operator with n = 4096, 
2048, 1024, 512 and 256 respectively. Figure 3.8 (MLE Segmented images of original 
and corresponding images approximated with De La Vallee-Poussin operator) 
contains MLE-based segmented images associated with the images in Fig. 3.7. The 
number of iterations taken for original image is 59400 and for n = 4096, 2048, 1024, 
512, 256 are 10850, 11850, 7500, 9450, 5400 respectively. 

The percentage LI and L2 percentage errors in the approximation for De La 
Vall6e-Poussin Kernel are summarized in Table 3.3. 

n 4096 2048 1024 512 256 

% Relative LI Error 9.52296 12.16614 14.63782 16.38696 17.40005 

% Relative L2 Error 9.62524 12.87752 15.87117 17.92870 19.08679 

Table 3.3: % Relative LI and L2 Errors in the Approximation of the Original Cross-Section with De 

La Vallde-Poussin Kernel 


In case of De La Vallde-Poussin operator for n = 128, 64, 32, 16 and 8, the 
smoothing increases which leads to a segmentation of interior brain region and its 
outer cranial part as two major segments. Interior part of the brain mainly represents 
two regions corresponding white matter and a combination of CSF and gray matter. 
With this operator also, the region represented by edema splits into two and then three 
parts as smoothing increases from n = 4096 to n = 256 and the boundary surrounding 
the edema thickens as smoothing increases. 



Original 


n=4096 



n=512 n=256 

Figure 3.7: Original T 2 weighted and Corresponding Images Approximated with De La Vall6e-Poussin Kernel 








n=512 


n=256 


Figure 3.8: MLE Segmented Images of Original and Corresponding Images in Figure 3.7 








CHAPTER IV : SEGMENTATION METHODS 

based on maximum likelihood estimation 


In this chapter, segmentation methods based on maximum likelihood 
estimation (MLE) have been studied. These methods have already been used in 
Chapters I, II and III in the numerical simulations. We have observed elsewhere that in 
MR images the distributions of pixel intensities of CSF, and normal tissues like white 
and gray matter, as well as different pathologies strongly seem to follow Gaussian 
distributions (R.K.S. Rathore, S. Datta, R.K. Gupta, S.B. Rao and R. Verma [2001]). 

Hence in the following treatment, the problem of segmentation has been 
modeled by considering the pixel intensities as random variables and the likelihood 
function defined as the joint probability density /(xr, 0, {j., a) of these variables, which 
is a function of the vector parameters 0 = (0i, 0^, ... , 0^)', \i. = {l^\, ixi, , pQ' and 
a = (cTi, 02 , • • ■ , Om)'. Considering the normal distribution of each tissue class within a 
pixel, the density function of the pixel intensity x is given as: 

, w dj {x-ixA^ 

J{x;9,pi,<j)= “Wexp^- - ^ >. 

j=\o 2a j 


4.1. INTRODUCTION 

The general method of maximum likelihood estimation is as follows: Let X 
denote the realized value of a set of observations and P(X,9) denote the joint density, 
where 0 = (0i, Oq), the vector of parameters belongs to a set 0 c Rq, where Rq is 
the q-dimensional real space. The likelihood of 0 given the observations is defined to 
be a function of 0: 


L(0|X) oc P(X, 0). 



The principle of maximum likelihood consists of accepting 0* = (©j", , 0 ^*) 

as the estimate of 0 , where 


L(0*|X) = supeg© L(0|X). 

In practice it is convenient to work with 1(0|X) = Log L(0|X), in which case 0* 
satisfies the equation 


l(0*|X) = supes 0 l( 0 |X). 


If the supremum is attained at an interior point of 0 and 1(0|X) is a 
differentiable function of 0 , then the partial derivatives vanish at that point, so that 0 * 
is a solution of the equations 


3/(0 J X) 
do I 


= 0 , 


/ = !, 


These equations are called maximum likelihood equations and any solution of 
them as a maximum likelihood estimate (C.R. Rao [1989]). 


Suppose M pixels on an image arc distributed into m segments with Oj as the 
;th population size ratio corresponding to the tissue class j (so that = !)• Then, 

the likelihood function for the totality of the pixel intensities, modeled as the mixture 
of the m number of normal populations is given by, 

fix-,0, n,a) = CTj)}, 

/=! 

where fj.j is the mean, Uj the standard deviation, and 6j is the size ratio of the jth 
population, x denotes the pixel intensity and the Gaussian density function is given by 


/V(x;/t,a) 



This single image case of the present model is discussed in the next section 
4.2, while the generalization of the same to the multispectral case of a set of images 
for the same cross-section is considered in section 4.3. The numerical and clinical 
details of the corresponding results in the test cases are presented in the later sections 
of the chapter. 


4.2. SINGLE CONTRAST SEGMENTATION METHOD 


In this section, we consider a single contrast normal distribution as follows: 

1=1 

The objective function log L is maximized with respect to = 1- In this 

case, the maximum likelihood function is obtained via calculating the estimators for 
which the likelihood function is maximum by solving the normal equations dUd^k = 
0, dL/dUk = 0 and dUdOk = 0, for /: = 1, ..., m, or by solving the normal equations d 
logUdiXk = 0, <? \ogUd(Jk = 0 and ^ logI/^0* = 0, for ^ = 1, ..., m, as both L(0) and log 
L{&) have maxima at the same values of Aik, oic and ic = 1, ..., m. Here /tk, and 6^ 
represent the corresponding mean, standard deviation and the population size ratio of 
the ^h segment present on an image. Let M denote the total number of pixels on an 
image. The following numerical iterative scheme (R.K.S. Rathore, S. Datta, R.K. 
Gupta, S.B. Rao and R. Verma [2001] ) is easily seen to be consistent with the 
maximum likelihood estimation (MLE) normal equations: 


Ml"*” = . 

(4.2.1) 

a<-‘> =/(l/£i.)f " , 

(4.2.2) 

e<"'' = (1/M)2" >iVU,;/,<”>,a<">)/S,. , 

(4.2.3) 

where 



(4.2.4) 

and 



(4.2.5) 


Iterative computation of means, standard deviations and population size ratios 
and pixel classification in the single image case, using above equations (4.2. 1-5) 
could be summarized as the following algorithm: 


Algorithm 4.2.1. 


Step 1: Given that m segments are present on an image, start with some initial 

estimates /^k*' \ ^ and \ k = 1, m, of their respective means, standard 

deviations and the population size ratios. 

Step 2: Update means [Ak^"\ standard deviations and population size ratios 

using equations 4.2.1-4.2.5, till ^ A: = 1,..., m, for some 

tolerance constant A > 0. (Note that k stands for a segment number and n for an 
iteration count.) 

Step 3: Calculate the probability function using the normal density function with 
respect to each of the segments for each pixel and classify the pixel to that segment 
for which the probability function is maximum. 


4.3. MULTI-RESOLUTION SEGMENTATION METHOD 

The contrast between the normal tissues like white matter, gray matter and the 
cerebrospinal fluid (CSF) is more on T 2 and Ti weighted images as compared to that 
on the proton density (PD) weighted images. T 2 weighted images show better contrast 
between white and gray matter, while Ti weighted images give a better contrast of 
CSF with that of the brain parenchyma (white matter and gray matter). Experience 
shows that it is not possible to get satisfactory segmentations using simple 
thresholding algorithms or similar such algorithms based on intensity only. This is 
because, in addition to the images not having a consistently good intra-tissue contrast 
across the entire boundaries, the presence of inhomogeneities within different tissue 
types results in the overlaps of different regions at their boundaries, as well as in the 
appearance of perforations in their interiors. The variability of the contrast in different 
tissue-types obtained by employing different sequences of magnetic resonance images 
suggests using multivariate analysis for the segmentation of a region of interest in a 
cross-section. 



Consider pixel intensity as a multivariate possessing different correlations 
between the pairs of given weighted images. The above methodology is not applicable 
in the present case and therefore the following iterative scheme is applicable to only 
cases comprising of two or more images. 

The derivation of the iterative scheme in this section is based on a multivariate 
normal distribution modeling of the mixutre of the tissue populations. From the 
definition of likelihood function, the model could be defined as: 

/=! 

where 

^ j^2'e^p|“ ~^j) 

The main aim is to estimate the mean vectors, covariance matrices and the 
population size ratio vector 9 = (du 62, , dm) which maximize log L subject to the 

condition = 1. This method is applicable to n number of images, where n is 

greater than one and is finite. In this model, x,- is a n x 1 vector for each pixel (i = 1, 
..., A/) on a cross-section, pi,, is an n x 1 vector for each jf = 1, ..., m, where m is the 

pre-assumed number of segments present on the cross-section and M is total number 
of pixels on it. Sy is an n x n positive semi-defmite covariance matrix for each j = I, 
..., m. We have considered the covariance matrices to be strictly positive definite. In 
this case also, the maximum likelihood estimators are obtained by solving the normal 
equations so obtained by using partial differentiations with respect to the estimators. 

The iteratively estimated mean vector multivariate distribution may 

be obtained by solving system of linear equations given by 

(4.3.1) 


where the elements of and for r,s=l,2, ...,n are as follows: 



where ^ d/ and 5^„^''^are the elements of matrix 

Syfc • At = 1, 2, . . ., m,, the values of 6 are obtained as 

(4.3,2) 

Due to positive definiteness of covariance matrix it can be 

decomposed into the product of a lower triangular matrix with its transpose using 

Cholesky decomposition. To compute the elements of covariance matrix we 

first estimated the elements of the lower triangular matrix of its Cholesky 
decomposition. The ratios of estimated off-diagonal elements to the estimated 
diagonal elements of the lower triangular matrix are obtained from the n-1 systems of 
linear equations given by 

=D^^''\^ = l,...,n-l, (4.3.3) 

where, for r, p > q, the elements of and D^^'^^are computed as follows: 

and 

The estimated diagonal elements for g = 1, 2, ..., n of the lower triangular 
matrix are obtained as 

(p (y) /v"' p 

^ (V) ^ ‘j. ' 

(4.3.4) 


Computation of mean vectors, covariance matrices and the population size 
ratios and pixel classification in the multi -resolution case could be summarized as the 
following algorithm: 



Algorithm_4.3.1. 


Step 1: Given that m segments are present on an image, start with some initial 
estimates of their respective mean vectors, covariance matrices and the population 
size ratios. 

Step 2: Estimate mean vectors, covariance matrices, and the population size ratios 
using equations (4.3. 1-4.3 .4), till = for some 

tolerance constant B > 0. (Here k denotes the segment, and v the iteration number.) 

Step 3: Calculate the probability function using the multivariate normal density 
function with respect to each of the segments for each pixel and classify the pixel to 
that segment for which the probability function is maximum. 


4.4. OBSERVATIONS RELATED TO TEST IMAGE SETS I AND II 

The segmented cross-section obtained using Algorithm 4.2.1 on test image set 
I with A = 0.001 is shown in Figure 4.1. The segmented cross-section obtained using 
Algorithm 4.2.1 on T 2 -weighted image of set II with A = 0.001 is given in Figure 4.2, 
in which ventricles and edema both fall under a single distribution. 

The segmented images so obtained on application of Algorithm 4.3.1. to all 
the possible pairs of images of image set II, with B = 0.02, are shown in Figure 4.3. 
Ventricles and edema merge to give a single segment in segmented image obtained 
using PD and T 2 weighted images, while proper segmentation of white matter and a 
part of edema adjacent to ventricles are found in segmented image obtained using PD 
and Ti weighted images, which also shows two different segments in ventricles. Only 
one segment represents both edema and ventricles in segmented image obtained using 
T 2 and Ti weighted images, which also shows a part of white matter merging with 
gray matter. 



T 2 weighted Image 

Figure 4.1: Segmented Image Obtained using Algorithm 4.2.1 to Test Image Set I 



T 2 weighted Image 

Figure 4,2; Segmented Image Obtained using Algorithm 4.2.1 to T 2 weighted Image of 

Test Image Set II 


Figure 4.4 shows a segmented image obtained using Algorithm 4.3.1 on the 
test image set II, with B = 0.02, taking PD, Ta and Ti weighted images. It is observed 
that PD, Ta and Ti weighted images together produce good segmentation of different 
norma! tissues along with the pathology present on the images. The segments 




corresponding to CSF, white matter, gray matter and edema are well represented on 
this segmented image. 



PD-Ta weighted Images 



PD-Ti weighted Images Tj-Ti weighted Images 

Figure 4.3: Segmented Images Obtained using Algorithm 4.3.1 to each pair of Test Image Set n 


Also the segmentation is done using three weighted images for consecutive six 
cross-sections obtained from a single patient and the results are shown in Figure 4.5. 
All the results are obtained with B = 0.02, and the segmented images so obtained are 




satisfactory. The number of iterations taken for cross-sections 1, 2, 3, 4, 5 and 6 are 
47 , 45, 34, 39, 50 and 59 respectively. 

With the variation in the convergence criterion, i.e., varying the value of B 
from 0.05 to 0.30 with an increment of 0.05, the results so obtained are presented in 
Figure 4.6. The number of iterations taken in cases of B = 0.05, 0.10, 0.15, 0.20, 
0.025 and 0.030 are 32, 14, 14, 14, 13 and 13, respectively. The sizes of major 
segments like white matter, edema, ventricles, etc. are almost similar in all these 
cases, but the graininess within the segments increases with the increase in the value 
ofB. 



Figure 4.4: Segmented Image Obtained using Algorithm 4.3.1 for PD, T 2 and Ti weighted Images of 

Test Image Set II 



Cross-section 1 


Cross-section 2 



Cross-section 3 



Cross-section 5 




Figure 4.5: Segmented Images of Six Cross-sections in Case of Edema 







B=0.30 


B=0.25 

Figure 4.6: Segmented Images of Test Image Set II in Case of Edema with Different Values of B 








B=0.30 




B=0.25 

Images of Test Image SetH in Case of Edema with Different Values of B 
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4.5. OBSERVATIONS RELATED TO TEST IMAGE SET III 

In test image set III, gliosis seems to be hyperintense on MT SE Ti-weighted 
image, and it is little visible on PD-weighted image. The contrast of either lesion or 
gliosis is not good on any of these weighted images. So, when we apply a 
segmentation algorithm based on MLE in case of a single image taking either PD, or 
MT SE T 1 weighted image, it gives single segment for gray matter and gliosis, where 
initial number of segments was fifteen. However, even with the increase in number of 
segments to more than twenty in the segmentation procedure, gliosis keeps on 
merging with the gray matter. This is because of the variation within the cranial part 
of the brain (consisting of skin, bone, muscle and fat) which is much more than that 
present between gray matter and gliosis, resulting in the further divisions of the 
segments within the cranial parts. 

Using multi-resolution Algorithm 4.3.1 to all possible pairs of PD, T 2 and Ti 
weighted images with B = 0.02, the results obtained are shown in Figure 4.7. In these 
cases, lesion and the perilesional gliosis are well separated in PD-T 2 pair, whereas 
white matter is segmented satisfactorily in PD-Ti pair. The result obtained applying 
multi-resolution Algorithm 4.3.1 with B = 0.02, taking PD, T 2 and Ti weighted 
images of test image set III, is also shown in Figure 4.7. The segmented image clearly 
shows the separate segments corresponding to lesion, perilesional gliosis and white 
matter. This has happened since the contrast between the lesion and the perilesional 
gliosis compared with the gray matter is high on T 2 weighted image due to T 2 
abnormality of pathology. So, when we combine PD, Ti, and T 2 weighted images for 
segmentation, it is observed that all major brain components like gray matter, white 
matter, CSF and abnormalities (lesion and gliosis) have distinct normal distributions. 
The classifications of different tissues were based on Bayes' decision theory (R.O. 
Duda and P.E. Hart [1973]), according to which the pixel belongs to that segment for 
which its probability is maximum among all the segments. 

The regions of lesion and ventricles are represented by same gray level on 
segmented image as both of them together form a single distribution. So, the 
distinction of the lesion from ventricles is difficult on the basis of gray level and hence 


PD-T 2 weighted Images 



PD-T 2 -T 1 weighted Images 


Figure 4.7; Segmented Images Obtained using Algorithm 4.3.1 to Test Image Set m 
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extraction of lesion has been done using seed - growing algorithm, with four or eight 
connectivity of pixels, having seed pixel from the region of the lesion. Similarly, also 
the white matter and the gliosis have been extracted out using seed growing algorithm 
with eight connectivity neighborhood. The extracted parts of the white matter, the 
lesion and the perilesional gliosis are shown in Figure 4.8. 


A class of different tissues may give one segment as in the case of lesion and 
cerebrospinal fluid present in ventricles, similarly one tissue may not make a single 
connected region like gray matter which comprises of several segments on an image. Two 



Figure 4.8: Extracted White Matter, Lesion and Perilesional Gliosis 





or more tissues falling under a single distribution can be extracted and the structure of 
gray matter comprising of different components can be obtained on the basis of 
anatomical knowledge of brain and the connectivity of the region. 


Graphs 4.1, 4.2 and 4.3 show the histogram plots of white matter, lesion and 
perilesional gliosis on PD, T 2 , Ti and MT SE Ti weighted images. Distribution of 
white matter shows normal on all the weighted images, while that of lesion is normal 
on T 2 weighted image and close to normal on other weighted images. Distribution of 
perilesional gliosis is close to normal on PD and MT SE Ti weighted images and two 
normal distributions are visible within gliosis on T 2 and Ti weighted images. Gliosis 
shows two peaks of normal distributions on T 2 weighted image, as it is represented as 
a T 2 abnormality on T 2 weighted image. Similarly gliosis is not separately visible on 
Ti weighted image as a separate entity. 



Graph 4,1: Histogram Plots of White Matter in PD, T 2 , Ti and MT SE Tj weighted Images 
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Graph 4.2: Histogram Plots of Lesion in PD, T 2 , Ti and MT SE Ti weighted Images 
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Graph 4.3: Histogram Plots of Perilesional Gliosis in PD, T 2 , T, and MT SE Tj weighted Images 


Means, variances and the correlation coefficients between the pairs of PD, Ta 
snd Ti weighted images in cases of lesion, gliosis, and their contra-lateral regions 
have been obtained after segmentation and the results are shown in Tables 4.1 and 
4.2. Results corresponding to white matter are separately given in Table 4.3. All the 
quantities are calculated on values normalized to 0-255. 


Image 

Lesion 

Contra-lateral 

PD weighted image (|x ± a) 

131.5617.54 

145.53 ± 6.77 

T 2 weighted image (p. ± a) 

192.03 ± 13.28 

127.63121.77 

T 1 weighted image (p ± a) 

68.53 ± 6.89 

106.76 1 12.57 

MT SE Ti weighted image (p ± < 7 ) 

64.85 ± 5.35 

75.1014.50 

Correlation coefficients (PD, T 2 ) 

0.068290 

-0.033513 

Correlation coefficients (PD, Ti) 

0.859128 

0.251560 

Correlation coefficients (Ti, T 2 ) 

0.148065 

-0.841327 

No. of pixels 

991 

991 


Table 4.1: Means, Variances, Correlation Coefficients of Lesion and its Contra-lateral part 


Image 

Gliosis 

Contra-lateral 

PD weighted image (p 1 a) 

153.7515.031 

141.4115.07 

T 2 weighted image (p 1 0 ) 

187.39 1 23.93 

120.82 1 14.32 

T 1 weighted image (p 1 a) 

102.54110.91 

107.06 1 10.49 

MT SE Ti weighted image (p 1 a) 

81.301 4.47 

74.47 1 2.95 

Correlation coefficients (PD, T 2 ) 

0.028229 

0.315359 

Correlation coefficients (PD, TO 

0.413414 

-0.073068 

Correlation coefficients (Ti, T 2 ) 

-0.790365 

-0.776473 

No. of pixels 

438 

438 


Table 4.2: Means, Variances, Correlation Coefficients of Perilesional Gliosis and its Contra-lateral part 


Image 

White Matter 

PD weighted image (p ± a) 

136.80 ±4.09 

T 2 weighted image (p ± a) 

111.94 ±7.45 

Ti weighted image (p ± a) 

112.01 ±3.26 

MT SE Ti weighted image (p ± a) 

73.45 ±2.52 

Correlation coefficients (PD, T 2 ) 

0.259644 

Correlation coefficients (PD, Ti) 

0.397903 

Correlation coefficients (Ti, T 2 ) 

0.172379 

No. of pixels 

7557 


Table 4.3: Means, Variances, Correlation Coefficients of White Matter 


From the computed values for weighted images in contra-lateral regions 
corresponding to lesion and gliosis, it is clear that they correspond to similar tissues in 
that. The Ta weighted values of lesion and gliosis are falling in one range, which is 
clear from the image also. Comparing the values in contra-lateral parts of lesion and 
gliosis and white matter, it becomes clear that the part of lesion along with gliosis 
falls within white matter. 

Segmentation using PD, Ta and Ti weighted images were obtained using 
Algorithm 4.3.1., with B = 0.02, for consecutive five cross-sections obtained from 
single patient, and the results are shown in Figure 4.9. The number of iterations taken 
for segmentation of cross-sections 1, 2, 3, 4 and 5 are 51, 59, 36, 45 and 43, 
respectively. The results corresponding to the variation in convergence critenon 
(varying the value of B from 0.05 to 0.30 with an increment of 0.05) are shown in 
Figure 4.10. The number of iterations taken in case of B = 0.05, 0.10, 0.015, 0.20, 
0.25 and 0.30 are 38, 17, 16, 16, 15 and 12, respectively. The sizes of major segments 
like white matter, ventricles, etc., are almost similar in all these cases, but the 
graininess within the segments increases with the increase in the values of B. Fo 
values of B from 0.05 and above, the parts of both lesion and penlesional gliosis have 

been merged with gray matter. 


Cross-section 1 


Cross-section 2 



Cross-section 3 


Cross-section 4 







B=0.25 


B=0.30 


Figure 4.10: Segmented Images of Test Image Set BI in case of Perilesional Gliosis with Different 

Values of B 
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In all our experiments, we have considered MR images of size 256x256. It is 
observed that the images reconstructed from 192x256 Fourier data are artificially 
extended to 256x256 grid by adding two strips of uniform distribution of zero values 
on both sides of the image. In such situations, Algorithms 4.2.1. and 4.3.1. were 
applied after the removal of these uniform distributions as both the algorithms would 
fail to give segments corresponding to uniform distributions due to zero variance in 
these distributions. 

Segmentation results into two or more number of segments due to the presence 
of noise in the outside region of the brain (mostly air) on the images, while it seems to 
possess a single segment. This can be suppressed by thresholding the image and then 
removing it after getting a uniform distribution outside, but it alters other segments in 
the cranial part also, as the intensity of bone matches with that of exterior part on all 
the weighted images. 

Sometimes segmentation leads to misclassification of certain pixels due to the 
inhomogeneity present within the tissues. This can be reduced by using a proper 
smoothing filter, which smoothens the region without disturbing the major edges 
before the application of segmentation algorithms or by considering further accuracy 
in the convergence criteria. 

Our further experiments show that if one removes the cranial part of brain, 
then one would get a better segmentation as the inhomogeneity within each tissue in 
the interior part of the brain is less as compared to that of the cranial part. Also, one 
can study the cranial part of the brain separately on the basis of this algorithm. 


4.6. BOUNDARY TRACING USING CUBIC SPLINES 


In most of the magnetic resonance images, the boundary of any tissue or 
pathology seems to be blurred due to partial volume effect. The method proposed here 
is helpful in tracing the boundary of any closed region comprising of pathology or any 
normal tissue. The boundary of a region on an image is calculated from the minimum 
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energy convex combination of the original image, its histogram and the Sobel edge 
map of it. This combination sharpens edge so as to enable an easy determination of 
knot points on it for the construction of the cubic spline (S. Datta, R.K. Gupta, S.B. 
Rao, V.S.N. Kaliprasad and R.K.S. Rathore [2000]). The 2-D cubic spline 
interpolationg curve zit) = (xit),y(t))' (J.H. Ahlberg, E.N. Nilson and J.L. Walsh 
[1967]) is obtained as 
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A binary image is obtained after tracing the boundary using cubic spline, 
which is further segmented into different components on the basis of a seed growing 
algorithm. The area inside the spline can be calculated by conventional method by 
integrating cubic spline and also by counting the total number of pixels comprising of 
particular tissue. 


The area inside the region is calculated by integrating the spline as: 


+ S",Kl/2)[(->:j -IXj’j-i -,,))) 

where Vyj — yj y ~ Xj — . (^,77) is any point inside the region. 




Original 


Boundaries of Tuberculoma Superimposed on 
Original Image 


Figure 4.11: Original Image and Boundaries of Tuberculoma Obtained using Cubic Splines 

Superimposed on Original Image 


Area under 

# Area 

J Area 

Outer curve 

1101 

1101.863 

Inner curve 

697 

697.28418 


Table 4.4: Total Area under Outer and Inner Rims of Tuberculoma 


#Area by calculating the average of no. of pixels inside spline and no. of pixels inside 
1 1^^ onH fArpfl calculated bv integrating spline curve. 
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This method is applied on an post contrast MT SE T i image which is obtained 
after injecting Gd-DTPA intravenously in the dose of 0.1 mmol/kg body weight. In 
this case the lesion is characterized by intracranial tuberculoma surrounded by edema. 
The inner and outer rims of tuberculoma are obtained using the proposed method. The, 
areas calculated inside the rim by integrating the spline and also by discretely 
calculating the number of pixels within region are comparable. The result is 
summarized in Table 4.4 and the original image and the image having boundaries 
superimposed on it are shown in Figure 4.11. The inner and outer boundaries obtained 
here are very much close to the boundaries, which are seen on original image. 

The method proposed here is useful in obtaining boundary of a lesion with 
discontinuous edges. This can be implemented for a boundary of any shape. Instead of 
storing all the boundary points, one can trace a complete boundary by storing knot 
points on the boundary. 


4.7. CONCLUSIONS 

It is not always that one has to use all of the three or more weighted images 
together to get a proper segmentation. This mainly depends on the region of interest 
and its contrast with respect to neighboring regions in each of these images, which 
suggests to select the number of images and the type of images to be used to get a 
proper segmentation. As we have discussed in section 4.4, to get the edema extracted 
from an image in case of image set I, it is sufficient to use only T 2 weighted image, 
while all of the PD, T 2 and Ti weighted images of test image set II were used to 
segment the edema, irrespectively of that both the image sets are of same patient with 
different cross-sections. In second case, even two images are not good enough to 
segment the region of interest. 

While in the case of post-traumatic epilepsy with perilesional gliosis, if one 
wishes to extract only gliosis part from an image, one can use only PD and T 2 
weighted images and PD and Ti weighted images are good enough for the 
segmentation of white matter. But, if we use PD, T 2 and Ti weighted images for 
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segmentation, then extraction of gliosis, white matter, lesion, etc., are possible after 
application of multi -resolution segmentation algorithm. 

The theoretical convergence of the MLE based iterative schemes is not known 
in the several cases that are possible. In all our experiments, nevertheless, we have 
observed the numerical convergence, except for the earlier mentioned cases of zero 
variances and with zero population size ratios, in 192x256 cases artificially extended 
to 256x256 grid. 
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